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ABSTRACT 


Three-dimensional (3D) woven composites have demonstrated multi-directional properties 
and improved transverse strength, impact resistance and shear characteristics. The 
objective of this research is to develop a new model for predicting the elastic constants, 
hygrothermal effects, thermomechanicai response, and stress limits of 3D woven 
composites; and to develop a computational tool to facilitate the evaluation of 3D woven 
composite structures with regard to damage tolerance and durability. Fiber orientations of 
weave and braid patterns are defined with reference to composite structural coordinates. 
Orthotropic ply properties and stress limits computed via micromechanics are transformed 
to composite structural coordinates and integrated to obtain the 3D properties. The 
various stages of degradation, from damage initiation to collapse of structures, in the 3D 
woven structures are simulated for the first time. Three dimensional woven composite 
specimens with various woven patterns under different loading conditions, such as tension, 
compression, bending, and shear are simulated in the validation process of this research. 
Damage initiation, growth, accumulation, and propagation to fracture are included in these 
simulations. 
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Chapter 1 


Introduction 

1.1 Three Dimensional Composites 

Three-dimensional composites are reinforced with three dimensional textile preforms 
which are fully integrated continuous-fiber assemblies with multi-axial in-plane and out-of- 
plane fiber orientations. These composites exhibit several distinct advantages which are 
not realized in traditional laminates. First, because of the out-of-plane orientation of some 
fibers, three-dimensional composites provide enhanced stiffness and strength in the 
thickness direction. Second, the fully integrated nature of fiber arrangement in three- 
dimensional preforms eliminates the inter-laminar surfaces characteristic of laminated 
composites. Third, the technology of textile preforming provides the unique opportunity 
of near-net-shape design and manufacturing of composite components and, hence, 
minimizes the need for cutting and joining the parts. The potential of reducing 
manufacturing costs for special applications is high. 

Three-dimensional textile preforms can be categorized according to their manufacturing 
technique. These include braiding, weaving, knitting and stitching. 
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Braiding preforms are formed with three basic techniques, namely two-step, four-step and 
multi-step braiding. In the case of two step braiding (Florentine 1992), the axial yams are 
stationary and the braider yams move among the axial yarns. In four-step braiding, all yam 
carriers change their positions in the braiding process and do not maintain a straight 
configuration. Multi-step braiding (Kostar and Chou 1994) is an extension to the four-step 
braiding. By allowing for both individual controls of a given track/column and the 
insertion of axial yarns, the range of attainable braid architecture is greatly broadened in 
multi-step braiding. 

In woven preforms, there are two major categories. The first is the angle-interlock multi- 
layer weaving technique which requires interlacing the yarns in three dimensions. The 
warp yarn in this three-dimensional construction penetrates several weft layers in the 
thickness direction, and therefore the preform structure is highly integrated. The second is 
the' orthogonal wovens, for which the yams assume three mutually perpendicular 
orientations in either a Cartesian coordinate system or a cylindrical coordinate system. The 
yams in the Cartesian weave are not wavy, and as a result, matrix rich regions often 
appear in these composites. 

The process of stitching is mainly based upon an existing technology for converting two- 
dimensional preforms to three-dimensional ones. This process is relatively simple. The 
basic needs include a sewing machine, needle and stitching thread. Major concerns of the 
stitching operation include depth of penetration of the stitching yams and, hence, the 
thickness of two-dimensional preforms that can be stitch-bonded, as well as the degree of 
sacrifice of the in-plane properties due to the damage to in-plane yarns. 
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Three-dimensional knitted fabrics can be produced by either a weft knitting or warp 
knitting process. The technique of knitting is particularly desirable for producing preforms 
with complex shapes because the variability of the geometric forms is almost unlimited. 
The large extensibility and conformability of the preforms enable them to be designed and 
manufactured for reinforcing composites subject to complex loading conditions. The 
versatility of knitted preforms offers a new dimension in textile structural composite 
technology. 

Even though three-dimensional (3-D) composites have attracted much interest due to their 
unique mechanical properties, such as enhanced transverse moduli and strength, and 
improved shear resistance and impact damage tolerance, the actual use of the 3-D 
composite material in engineering design poses many problems. The main problem comes 
from the complex geometry of 3-D composites. The fiber geometry is so complex that the 
geometric modeling itself is very difficult, much less accurate stress distributions. For 
example, in plain weave textile composites, there are many fiber tows (warp and fill) 
interlacing each other. There can be nesting of the fiber tows of one layer in adjacent 
layers. The existence of matrix pockets adds to the complexity of the geometry. In fact, 
many research papers have been devoted just to modeling the geometry of 3-D composites 
(Pierce 1987, Pastore and Ko 1990, Du and Chou 1991). 

The inherent geometric complexity of 3-D composites makes a detailed stress analysis 
very formidable. Most analytical and numerical techniques are merely used to predict the 
stiffness properties of 3-D composites. There are few models that have been developed for 
detailed stress analysis (Woo and Whitcomb 1994, Yoshino et al 1981, Kriz 1989, Pastore 
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et al 1993) and the strength prediction of textile composites. To date, there is no 
information in the existing literature on simulating the entire procedure of damage 
propagation of 3-D composites. Some of the earlier works for predicting elastic constants, 
damage propagation, and strength of woven and braided composites are highlighted 
below. 

1.2 Stiffness Properties of Three Dimensional Composites 


The stiffness averaging method which was developed by Kregers et al (1978, 1979) is 
widely used to predict the deformation characteristics of a composite with three- 
dimensional reinforcement from the known mechanical properties of its components. The 
basic idea behind stiffness averaging method is to treat the fibers and matrix as a set of 
composite rods having various spatial orientations. The local stiffness tensor for each of 
these rods is calculated and rotated in space to fit the global composite axes. The global 
stiff n ess tensors of all the composite rods are then superimposed with respect to their 
relative volume fraction to form the composite stiffness tensor. This approach is also 
called the Fabric Geometry Model (Pastore and Gowayed, 1994) or Orientation 
Averaging Method. The stiffness of the individual directions of reinforcement are averaged 
in accordance with the following expressions: 
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Where A/« m are the components of the stiffness tensor of the three-dimensionally 
reinforced composite: V. is the calculated volume of the i-th direction of reinforcement; N 
is the number of discrete directions of reinforcement, N > 1. 

Ko (1986) presented a geometric model for three-dimensional braid composite using the 
concept of average cosine to evaluate the tensile strength and modulus of a three- 
dimensional braid composite. The three-dimensional braid composite was divided into a 
series of unit cells and the yarn segment was idealized as straight in a unit cell as shown in 



Figure 1.1 Idealized fiber yarn segment in unit cell of 3D braid composites 

To obtain an average representation of yam orientation, the average cosine of yarn angles 
was used: 

cos 6 = ND y / D f 
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Where N = number of yams in the fabric, D y = yarn linear density, Df = fabric linear 
density. The predicted composite tensile strength and modulus in the longitudinal direction 
based on the average angle of yam, in general, were within 20% of the experimental 
results. 

To facilitate the description of complex geometry of 3-D composites, the concept of unit 
cell structure has been used by many investigators. A unit cell is a representative volume 
element small enough to reflect the structural details, yet large enough to represent the 
overall response of the composite. When unit cells are repeated, they will reconstruct the 
entire structural geometry. 

Yang et al. (1986) presented a “Fiber Inclination Model” to predict the elastic properties 
of 3-D textile structural composites. The unit cell of a composite was treated as an 
assemblage of inclined unidirectional laminate. For a 3-D braided composite, the 
orientation of the yarns in the braided preform is controlled by the three orthogonal 
motions. Therefore, the resultant preform is a continuous interwoven structure composed 
of yams oriented in various directions. An idealized unit cell structure is constructed based 
upon the fiber bundles oriented in four body diagonal direction in a rectangular 
parallelepiped which is shown schematically in figure 1.2, the unit cell is considered as an 
assemblage of four inclined unidirectional laminae, and each unidirectional laminae is 
characterized by a unique fiber orientation and all the laminae have the same thickness and 
fiber volume fraction of each lamina is assumed to be the same as that of the composite. 
The laminate approximation of the unit cell structure is shown schematically in figure 1.3. 
The properties of each lamina can be obtained from the classical laminated plate theory, 
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the composite properties can be obtained from these lamina properties based on the 
stiffness averaging approach. 



Figure 1.2 The unit cell structure of a 3-D braided structural composites 
with yarns moving in three orthogonal directions 



Figure 1.3 The unit cell of the “Fiber Inclination Model” composed of four 
unidirectional laminae 
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Whitney et al. (1989) extended the Fiber Inclination Model to predict the in-plane elastic 
properties of 3-D angle-interlock textile structural composites. In the case of angle- 
interlock geometry, the unit cell is inherently more complex. Since there are no repeatable 
units in the thickness direction, the unit cell essentially occupied the entire preform 
thickness. Also, the yams can occupy any number of positions within the unit cell. Yams 
in warp and weft directions may occupy positions which are horizontal or inclined in the 
thickness direction. To account for varying bundle positions, the authors divided the unit 
cell into micro-cells that are repeated to reconstruct the entire unit cell, and the 
calculations are made on a generic unit cell with micro-cells. Crimping of fibers at the 
corners of the cell and the intersection of fibers at the unit cell center are ignored. 

Stiffness averaging method was also used to predict elastic constants of two-step braided 
composites by Byun et al. (1991). The architecture of this material was investigated by 
identifying the geometric and braiding process parameters which include the linear density ' 
ratio between axial and braider yarns, the aspect ratio of axial yam and the aspect ratio of 
braider yarn. The predicted results of axial tensile modulus shows reasonably good 
agreements with test results, however, shear moduli show a lack of agreement between 
prediction and measurement. 

All the aforementioned models are based on the concept of averaging stiffness. Ma et al. 
(1990) developed another methodology based on an energy approach to evaluate elastic 
stiffness of 3-D braided textile structural composites. Three types of elastic strain energies 
in the composite rods were taken into account. These included the strain energies due to 
bending, extension, and compression over the region of fiber contact. Closed form 
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expressions for axial elastic moduli and Poisson’s ratio were derived as functions of fiber 
volume fraction and fiber orientations. Numerical results show that the axial elastic moduli 
are sensitive to the geometrical braiding pattern. The moduli increase as the yam 
orientation angle decreases, that is, the fiber becomes more aligned with the tensile axis. 

Fiber bundle orientations are the essential geometrical properties for determining 
mechanical behavior. For the sake of simplicity, the geometry of 3-D composites have 
been represented by idealization rather than reality by most researchers. Most models have 
dealt with perfectly regular structures. The fibers are assumed to be straight inside a unit 
cell. However, real composites are highly irregular, and waviness of fibers is unavoidable. 
The angles between axial tows and braiders change in an irregular way. In order to take 
account of these irregularities, Cox and Dadkhah (1995) promoted a more practical 
approach to characterize the geometrical inconsistencies of typical triaxially braided 
composites. Out-of-plane misalignments were taken into account by introducing a 
cumulative probability distribution (CPD) for the misalignment angle £ . Typical CPDs can 

be fit quite well by a symmetric normal distribution F ^( %) with corresponding density 

function (£ ) = dF, / d £ given by 



1 



exp(-^ 2 /2cr^). 


The standard deviation <7^ can be determined experimentally. The Young’s modulus E x 
can be determined by 
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Where E x (£) is Young’s modulus for a unidirectional composite under a load oriented 
at angle £ to the fiber direction X. E x (q) can be expressed by 


1 cos 4 £ 




+ cos 2 £ sin 2 E, 




V G ,v ^ , 


+ 


sin 4 E, 


1.3 Stress and Strength Analysis of Three Dimensional 
Composites 

Strength of 3-D composites is more difficult to predict than stiffness. This section 
summarizes work done on predicting the strength, fracture, and damage tolerance 
behavior of 3-D composites. 

Ishikawa and Chou (1983) developed a 2-D fiber undulation model based on classical 
lamination theory which considered the undulation in the fill yarn direction but neglected 
the undulation in the warp yarns of woven composites. They modeled nonlinear shear 
response of both the fill yarns and the interstitial matrix along with the effects of 
transverse cracking of the warp yarns to predict the knee in the stress-strain response of 
woven composites. Kriz (1985) used a generalized plane strain finite element model which 
assumed a linear undulation path for the fill and warp yams to study the effect of 
transverse cracking on the stiffness and internal stresses of a glass/epoxy plain weave 
composite. 
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Stanton and Kipp (1985) developed a nonlinear constitutive model for plain weave 
carbon-carbon composites which accounted for the differences in tension and compression 
response, warp-fill crossover damage and multiaxial stress interactions using a simple 
interaction formula along with experimental stress-strain data for all six components under 
tension and compression loading. Jortner (1986) developed a 2-D mechanistic model 
which modeled the undulations of the fill yarns but neglected the undulations of warp 
yarns. He used a stress-averaging scheme which accounted for the nonlinear shear 
response of the fill yams and the transverse cracking of the warp yarns to analyze plain 
weave carbon-carbon composites. 

Ko and Pastore (1985) used a fabric geometry model to define the yam orientations in a 
three-dimensional braided composite. They used the yam orientations to first estimate the 
strength of the fabric preform and then compute composite strength using a simple rule of 
mixtures. Ko (1989) also used the fabric geometry model together with a maximum strain 
energy criterion to predict yam failures and ultimate strength of a 3-D braided composite. 

Dow and Ramnath (1987) modeled woven fabric composites using a simple geometry 
model that assumed a linear undulation path for the fill and warp yarns. They computed 
constituent fiber and matrix stresses from local yarn stresses which were calculated using 
an iso-strain assumption and predicted failure based on the average stresses in the fiber 
and the matrix along with a maximum stress criterion. Dow and Ramnath (1987) used the 
fabric geometry model together with a simple linear yam bending model and an iso-strain 
assumption to compute average fiber and matrix stresses which were used to predict local 
yam failure and strength for 2-D triaxial braided composites. 
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Naik (1994 & 1995) developed a micromechanics analysis technique for the prediction of 
failure initiation, damage progression, and strength of 2-D woven and braided composite 
materials. The yam architecture was discretely modeled using sinusoidal undulations at 
yarn crossovers. Overall thermal and mechanical properties were calculated based on an 
iso-strain assumption. This analysis technique included the effects of nonlinear shear 
response and nonlinear material response due to the accumulation of damage in the yams 
and the interstitial matrix and also the effects of yarn bending and the geometrically 
nonlinear effects of yarn straightening/wrinkling during axial tension/compression loading. 
Parametric studies were also performed on the woven and braided architectures to 
investigate the effects of yarn size, yarn spacing, yam crimp, braid angle, and overall fiber 
volume fraction on the strength properties of the textile composite. 

Three-dimensional finite element models (FEM) have also been used for the prediction of 
nonlinear material properties (Bhandarker et al 1991) and the modeling of damage 
(Blackketter 1993) of plain weave composites. The 3-D FEMs are highly computer 
intensive and also require considerable time and effort for model generation. Foye (1993) 
developed homogenized replacement finite elements and analyzed sub-cells within the 
repeating unit cell to overcome these limitations. However, he had to manually calculate 
the orientations of the different yam directions, yam interfacial planes and volume 
fractions, for each element in the 3-D model. FEM is, therefore, not well suited for 
performing parametric studies to investigate the sensitivity of composite strength 
properties to fabric architecture parameters. 
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Compared with the regular finite element analysis, a global/local finite element method 
which was developed by Sun and Mao (1988 & 1990) is more efficient to analyze textile 
composites. The global/local finite element method was based on the regular finite element 
method in conjunction with three basic steps, i.e. the global analysis, the local analysis and 
the refined global analysis. In the first two steps, a coarse finite element mesh was used to 
analyze the entire structure to obtain the nodal displacements which were subsequently 
used as displacement boundary conditions for local regions of interest. These local regions 
with the prescribed boundary conditions were then analyzed with refined meshes to obtain 
more accurate stresses. In the third step, a new global displacement distribution based on 
the results of the previous two steps was assumed for the analysis, from which much 
improved solutions for both stresses and displacements were produced. 

Woo and Whitcomb (1994 & 1996) utilized global/local finite elements to analyze textile 
composites. A relatively crude global mesh was used to obtain the overall response of the 
structure and refined local meshes were used in the region of interest where rapid stress 
change may occur. The homogenized engineering properties could be used in the crude 
global mesh away from the free boundary and transitional regions. In the local meshes, 
the details of the coarse microstructure of textile composites(e.g., the individual tows and 
matrix pockets) could be modeled discretely. In the transitional range of microstructure, 
however, discrete modeling may not be practical even with supercomputers due to the 
huge computer memory and CPU requirements. Use of homogenized material properties 
is also usually inappropriate. In this range, special macro finite elements can be used. The 
macroelements employ a number of subdomains or subelements to account for the 
microstructural details inside individual elements. Whitcomb et al. (1994) discussed two 
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types of macroelements. The elements described in these references are based on single or 
multiple assumed displacement fields. 

The global/local method with local refinement assumes that the region which requires 
further mesh refinement can be localized. If the solution behavior is complicated 
everywhere, the correction by the global/local iteration does not necessarily reflect the 
nature of the true solution accurately. For example, if the solution behavior is complicated 
everywhere and only some portion of the global domain is refined in the local mesh, the 
global/local iteration may do more harm than good in solving the problem. A solution to 
this problem is using an engineering global/local analysis. In this procedure, the initial 
coarse global solution is assumed to be close enough for the purpose at hand. That is, no 
global/local iteration is employed. The local problem is solved only once with boundary 
displacements from the coarse global solution. 
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Chapter 2 


Analytical Characterization 
of Three Dimensional 
Woven Composites 


Composite materials reinforced with woven fabric preforms are being considered for 
potential structural application in the aircraft and automotive industries. Fabric-reinforced 
textile composites potentially have better out-of-plane stiffness, strength and toughness 
properties than conventional laminates. However, the architecture of a textile composite is 
complex and, therefore, the parameters controlling its mechanical properties are 
numerous. Mechanical testing to characterize all the effects of textile architecture could be 
an economically unrealistic proposition. In order to facilitate the design of fabric- 
reinforced composite structures, it is crucial to develop experimentally verified analytical 
models for the prediction of both the uniaxial and multiaxial stiffness and strength 
properties. 
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2.1 Woven Patterns 


Xu et al. (1994) provided complete statements of the rather complex sequencing of 
through-the-thickness yarns. Figure 2. 1 shows the three typical types of weave in woven 
composites, namely, (a) layer to layer, (b) through the thickness angle interlock, (c) 
orthogonal interlock weaves. 

In woven composites, the stuffers and fillers alternate in layers through the thickness. The 
stuffers and fillers form a coarse 0°/90° array shown in figure 2.1 for most woven 
composites. Nevertheless, the proposed method allows the arbitrary orientation of stuffers 
and fillers in the X-Y plane. The through the thickness reinforcement, or warp weavers, 
may be oriented in any direction with reference to the 3D composite coordinate axes. 
Stitched composites may be modeled by weaver or stitch fibers that are oriented 
perpendicular to the X-Y plane and parallel to the Z-axis of the composite. 

2.2 Fiber Arrangement 

Complete understanding of arrangement of fibers, including stuffer, filler, and warp 
weaver fibers, is essential to predict the engineering properties accurately. Woven 
composites are typically composed of n s layers of stuffers and n s + 1 layers of fillers through 
the thickness as shown in Figure 2.2. 

In order to estimate the elastic properties of 3-D woven composites, an approach based on 
modified laminate theory and orientation average is used. First of all, composites are 
divided 
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into plies so that every ply contains either stuffer or filler fibers as shown in Figure 2.3. 
Since the weaver fibers go through the entire thickness, each ply must contain some 
weaver fibers for woven composites. 

For orientation averaging, a three-dimensional fiber arrangement can be considered as 
composed of two types of structures, the baseline structure formed by fibers in X-Y plane 
direction, such as stuffer fibers and filler fibers, and the interwoven or braided structure 
composed of warp weaver fibers which penetrate the thickness of composites. Therefore 
each ply contains two types of domains, the primary domains which consist of baseline 
fibers, and the weaver domains which consist of weaver fibers, the occupying volume 
fractions of these domains are denoted by f p and f w respectively. Each domain is 

characterized by an orientation along which the fibers within it are presumed to lie. All of 
the fibers in the primary domains are assumed to be parallel to X-Y plane and each ply can 
consist of the primary fibers in the same direction. Warp weaver fibers are always assumed 
to be piecewise straight and their orientations are defined by the directional angles to the 
X,Y,Z axes. 

Let [£'"’] denote the stiffness matrix for domain a. Based on orientation averaging 
method, the composite stiffness matrix [E c J of 3-D woven composite can be predicted by 
the following equation: 

W=I/j£“] (2-1) 

a 

Where [E un ] denotes [£ ( “ l l transformed from domain a material coordinate system into 
the composite coordinate system and f a denotes volume fraction of domain a. 
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i i Stuffer 
^^3 Warp weaver 


Figure 2. 1 Schematics of (a) layer-to-layer angle interlock, (b) through-the- thickness 
angle interlock, and (c) orthogonal interlock weaves. 
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Figure 2.2 Representative layer sequence of fillers and stuffers through the thickness, 
with the layer thickness tj and t s defined. For the case shown, n s = 2. 



Figure 2.3 Schematic of ply division in 3D woven composites 
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2.3 Three Dimensional Hygrothermoelastic Properties 


The first step in the computation of composite properties for any fiber configuration is the 
evaluation of the local orthotropic properties of a unidirectional composite ply. From 
composite mechanics (Jones, 1975), the compliance (inverse stiffness) matrix in material 
axes for each orthotropic ply written in terms of the engineering constants is 


p.r= 
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( 2 - 2 ) 


Total ply stiffness properties of a composite layer containing one primary in-plane fiber 
orientation and N w weaver orientations can be obtained from combination of the primary 
domain and weaver domain properties according to their volume fractions. Note that it is 
possible for a ply to have multiple weaver domains, but only one primary domain is 
considered. Therefore the stiffness matrix for the ith layer containing the primary and 
weaver domains may be written as: 


M = /,[*.- T[ £ 'l [*<']+ E /.,[<'] [ £ ."'] [*,''] 

2=1 


(2-3) 


where 
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£, is the three dimensional lamina stiffness matrix of ith ply. 


R) P is the two dimensional coordinate transformation matrix from primary domain 
material axes to composite structural axes. If the angle between material axes and 
structural axes is 0 , R? can be expressed by 


R; P = 


cos 2 # 
sin 2 # 
0 
0 
0 


sin 2 # 
cos 2 # 
0 
0 
0 


-cos# sin# cos# sin# 


0 0 0 

0 0 0 

1 0 0 

0 cos# -sin# 
0 sin# cos# 

0 0 0 


2 cos# sin# 

- 2cos#sin# 
0 
0 
0 

cos 2 # -sin 2 # 


/ is the primary domain volume fraction 

N w is the number of weaver domains, each weaver domain consists of weaver fibers in 
the same direction. If weavers contain fibers in different directions, they should be 
grouped into different weaver domains. 

is the three dimensional coordinate transformation matrix from jth weaver domain 
material axes to composite structural axes. Since directions of weavers are 
arbitrary, three dimensional coordinate transformation matrix is used here. If the 

direction cosines between the weaver local axes and structural axes are //, nii 

rij, R*' can be written as (Lekhnitskii 1977), 
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(n 2 l 2 +n y 

up 

lnhpx x 

ln 2 n x 

(j^h+Ph) 

[m 2 n ] +m x n 2 ) 

(np +n x i 


f w is the jth weaver domain volume fraction 

N w 

Note: 4 +X/h-,. =1 
j = i 


Ep is the stiffness matrix in the local coordinate system of jth weaver in ith ply. Before 
damage occurs, Ep should be same for each / since we assume every layer to have the 
same weaver arrangement. But after damage, every layer may have its own property 
degradation which could be different between layers, therefore degraded Ep could be 
different for different i. 

In order to differentiate three dimensional properties from two dimensional properties 
which are discussed in the next section, an overbar is added to three dimensional stiffness 
matrices and coordinate transformation matrices, such as [e,] and [/?,■]. 

The averaged properties of the entire composite laminate can be obtained by summing up 
all the layer properties and the properties due to interply distortion energy (Murthy and 
Chamis 1986). The interply distortion energy terms represent the stiffening effects of 
changes in the primary fiber orientation angles between adjacent layers. Therefore the 
composite laminate stiffness matrix is written as: 
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where 

E c is the three dimensional stiffness matrix of entire composite. 
t c is the thickness of the entire composite. 

N, is the number of plies of the entire composite. 

Z; is the distance from reference plane to the bottom of ith ply. 
H j is the interply layer distortion energy coefficient. 


Sj is generated from the following equation (Murthy and Chamis, 1986): 

A 2 - A 2 0 0 0 -AB 
-A 2 A 2 0 0 0 AB 

0 0 0 0 0 0 

(2-5) 

0 0 0 0 0 0 

0 0 0 0 0 0 

-AB AB 0 0 0 B 2 

A = sin 2 6 j+] - sin 2 

B - cos 2 6 j+] - cos 29 j 

Bj is the fiber angle of jth ply from material axes to composite structural axes. 

Thermal and hygral coefficients of expansion about the composite structural axes can be 
obtained from the following equations: 
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l C i=l 7=1 


( 2 - 6 ) 


l C 1=1 7=1 


( 2 - 7 ) 


Where {a c },{/?<-} are thermal and moisture coefficients of expansion about the 
composite structural axes. 


K } — {«C, , ’ «C 22 ’ «C„ ’ «C, ’ «C„ ’ «c 12 } 

{A; } = {Pq , ’ Pc 21 ’ Pc„ > Pc 2 , > Pq, > Pq 2 } 


( 2 - 8 ) 

( 2 - 9 ) 


2.4 Two Dimensional Hygrothermoelastic Properties 


In most cases computational evaluation of composite structures is carried out by using 
plate and shell finite elements. Therefore, the composite properties need be expressed as 
two-dimensional plate stiffness with membrane and bending degrees of freedom. The 
resultant forces and moments in terms of the middle surface extensional strain and 
curvatures can be expressed as: 



'[ a ] M1JM1JK}! 
Ib] mJImJ |{M r }j 


Kif 


( 2 - 10 ) 


where 
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{N] = {N x ,N y ,N xy } T 
{M} = {M x ,M y ,M xy } T 

{e} = {e x ,e,,e v } is the reference plane membrane strains. 

{k} - |?c v , K y , K xy } is the reference plane curvatures. 

Wr}. Kl is the resultant forces due to temperature and moisture change. 

{ M T }, { M M } is the resultant moments due to temperature and moisture change. 

For a three-dimensionally reinforced laminate containing layers of primary and woven fiber 
reinforcement domains, the resultant forces due to external load can be obtained from the 
following equation: 

Ki=ir fri* 

;=i i,_1 

= £fi [£,]«£}+ *{«■})& 

,=i 

= £ f' 4 [« ' f k Ikl + £ /.. ]]({e}+ z{k})cIz 

'■=> ‘ ,_1 L M 

(2-11) 

where 

[E; ] = /, k f k ]k ] • + £ r., k ] <2-12) 

7=1 
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t A ']=X( z . 
1=1 


1=1 


[®r'] 

j =l 

4[ R '] r [ £ / , ][«.'']+i/.,[ g «''] 


and E i is the two dimensional lamina stiffness matrix of ith ply. 


E? is the two dimensional lamina stiffness of primary domain of ith ply in local axes. 



1 
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v 2l 



0 


0 

0 
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K is the two dimensional coordinate transformation matrix of primary domain. 


cos 2 0 , 

sin 2 0. 

. —sin 20. 


( 

2 

sin 2 0, 

cos 2 0. 

sin 20 ; 

t 

t 

2 

- sin 20 ; 

sin 20, 

cos 20, 


E™' is the two dimensional stiffness of weaver domain in structural axes. It can be 
obtained by extracting relevant terms, which correspond to On, G 22 , O 12 , from the three 
dimensional stiffness matrix that can be generated by the following transformation: 

S'- 

Superimposing the stiffness contributions due to interply layer distortion energy: 
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(2-13) 


N,~ 1 


[4 = [A']+£ff,[sJ 

j= i 

-ik - v, )Lkn^'] [^]+£/.,[«.'']|+"f «/W 


N,-l 


[»] = [B'] + S<W,[S,] 

7=1 



(2-14) 


where z' is the distance from the interply layer to the composite reference plane (mid- 
thickness). 


Similarly, 


{ M J = X£ t W,}zdz 

= ]£ f [ £ « ]K e } + *M) z dz 

i=i z '-> 

=if 4 fr'ri £ '] [*']+£/., [«,“■' ] (H+iWii* 


=R'M tW 4 


where 


[°'H i t? - zt. f 4 k r k k /., [ e ? i 

i=l >1 


(2-15) 
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Superimposing the stiffness due to interply layer distortion: 


N ,- 1 


L 1=1 




J i= 1 


4kr[ £ /]kl+i/., [£,"'] +\%?H,[s l ] 

M z 7=i 


For the primary domain, the thermal forces are 
K' } = £ k ] AT {«,'}* = (z, - Z H ) AT k ]{«;} 


For the yt/z weaver domain, the thermal forces are 
K' } = £ [K 1 } Ar fr' }* = (z i -z.-_ 1 ) Ar [Ep ] {ap } 


( 2 - 16 ) 


( 2 - 17 ) 


( 2 - 18 ) 


Transforming Nj and Np to global axes, and then averaging them in terms of the 

weight of their volume fraction to get resultant force for ith ply in global coordinate 
system due to temperature change: 


W}= /, k IK } + i k ]k? } 

i = 1 

=A7;( z ,-z w )|/,ki£.'K“'} + i/.,k , ''K'K'}] 


7=1 


( 2 - 19 ) 

Summing up the Nj for all phes to get the resultant thermal force for entire composite 
structure: 
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Rl = £{*>■,} 


/= 1 
N. 


/= 1 


at;. (z, - z H K/, ] [# IK } + i /„, ] [e? ]{ ccp } 


Similarly, {a m } , { M r } and {m m } can be generated by: 


( 2 - 20 ) 


Ri=i 


i=l 


AM, (z, - z H )| 4 M*,' ]{A P } + 1 U, [*? ] [e? ]{pp } 


( 2 - 21 ) 




2f. 


at; (zf - zl)\f p [R!’}[E!’]{a!’} + 


( 2 - 22 ) 


{"-Kx 

z i=i 


AM,, (zf - zU )|/, [R?][E' ]{A' } + £/„, [<' ][C J ]{/T } 


(2-23) 


2.5 Stresses in the primary domain and in the weaver 
domain 

In order to perform damage propagation simulation, one needs to evaluate the stresses in 
stuffer, filler, and weaver fibers and matrices. For most composite structures, the first step 
for stress analysis is the assembly of a finite element model using the stiffness properties 
given in the previous section. After finite element analysis, the resultant forces and 
moments are obtained for each node. From equation (2-10), we can write: 
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(2-24) 



[A] [*]' 

[B] [D] 


The strains and stresses in global axes for a specific ply — the ith ply can be obtained from 
the following equations: 


{e,} = {e}-zW 

= [A]" {{/V } + K } + {/V M } - [B]{k}} -z{k} 


{o', } = [E, ] ({e, } - A 7) {a, }- AM, {fi }) 


(2-25) 

(2-26) 


where 

tel = 1^.4'. 4'1 

(o', } = {o-;; 1 . o-;; 1 , a-;; 1 } 


the other stress components ,o [ ^ can be obtained from the transverse shear 

resultants Qn, Qn and the transverse pressures P u , P, by the following equation 
(Vladimir, 1975): 


Ict'V 


33 





16 z, 

Iff 




(2-27) 
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Stresses in the primary and weaver domains are computed from the corresponding strains 
using the appropriate stress-strain relations. Since the primary domain is in the same plane 
as the ply, it has the same strain components as the ply, that is 

} = [#]{*,} (2-28) 

and the primary domain stresses are computed from: 

k } = k } W; ] [A]- 1 (m + k } + {n m } - [b]{k}) 

- [£'](A7;.{ar } + A M ,{&'}+ z[r, p ]{k}) 

(2-29) 

Since the weaver is not in the same plane as the ply, one needs to transform the three 
dimensional strain components of the ply to the weaver axes to obtain the weaver strains. 
Weaver domain stresses can be computed from the weaver strains. 

The three dimensional stress-strain-temperature-moisture relationship can be expressed as: 
kJ = K](k - A7 ;{«;} - AM, {A}) (2-30) 

where 


kl = k; 1 ’ 
kkk' 


'-'22 ? 1^33 ’ '-'23 ’ 


pin p('> p<i> 

, c-22 1 *-33 ’ *-23 9 


(i) 


OV3> 


c (n 

-13 » 



From (2-25) and (2-27), we computed £,‘ 2 \ <7,^ 1 , cr , cr^ V • Now, in terms of (2- 

30), we may find the other three strain components e',,' . Transferring the strains 

{fi, } to weaver axes we obtain the weaver strains in local material coordinates: 
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(2-31) 


Then the weaver stresses can be expressed as: 

k’ } = k ] (k } - A T, } - am, {/f }) (2-32) 

2.6 The Influence of Fiber Waviness 

In woven 3-D composites, nominally straight in-plane yarns, stuffers and fillers, are often 
much more wavy than those in conventional laminates because of the existence of weaver 
fibers. The most important effect of tow waviness on elastic properties is to reduce the 
axial stiffness of a tow. Since the microstructure of a woven composite is highly complex 
and the waviness of fibers is random rather than possessing a single characteristic 
wavelength or amplitude, it is not currently practical to model the actual tow geometry. 
Idealization is mandatory. For simplicity, a composite unit cell with wavy fibers shown in 
figure 3.4 is analyzed in this section. 

To simulate the fiber configuration, the following form of waviness is assumed, 

y-~^-(lx-x 2 ) (2-33) 

Where l is the half-wavelength and C is the amplitude of waviness. Therefore, the length 
of yam S(l,C) can be expressed by 
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s(l,c)=j‘*Jl+y ' 2 dx 
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(2-34) 


/ 2 


l 



Figure 2.4 Wavy tow model for analyzing the influence of fiber waviness 
on longitudinal stiffness of unidirectional composites 

Supposing that the composite unit sustains axial tensile load Gi as shown in figure 2.4, the 
extension of the unit in the X direction is caused by two parts: one is due to fiber axial 
extension under axial tensile load, and the second is due to the reduction of fiber waviness, 
i.e. the fiber will become more straight under the axial tensile load, and therefore, 
additional extension will occur to the composite unit. Since the first component has been 
widely studied in the literature, this section is focused on the second component— the 
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extension due to the reduction of fiber waviness. Fiber axial extension is not taken into 
account when considering the second component contributing to composite extension. 
Therefore, it is reasonable to assume that the length of fiber yarn remains unchanged, i.e. 
5(7, C) is constant when analyzing the second component of extension. When the fiber 
waviness reduces, the amplitude C decreases and the wavelength increases. 

Since 5(7, C) is constant, one may write 
5(/,C)=|y<5/+||sC=0 

Therefore the change in composite length becomes: 

SI=-[^I^\SC=-F{UC)SC (2-35) 

^ u C f O t j 

Where 


F(l,C)= 



(2-36) 


Assuming the original length of the composite unit is l 0 , after the axial load is applied, the 
length increases to /, then extension of the unit can be expressed by 

(-'o=|H+<5/ (2-37) 

E \\ 
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Where £j° is the composite longitudinal modulus based on orientation averaging method. 

— £-/ 0 is the extension due to fiber axial extension, 81 is the extension due to reduction of 
E u 

fiber waviness. 

When the composite unit cell is subjected to tensile load T c -(J l A c , the fiber bundle is 
subjected a tensile force T f -a f A f . Since the fiber bundle is not straight, there must be 

balancing reaction q as shown in figure 2.4 to counteract the vertical component of T f . 
According to equilibrium equation, q can be written as 

d 2 y 8C 

q--a f Af—^--jY' < jfAf (2-38) 


Where o f is the longitudinal stress in the fiber and A f is the cross section area of the 


fiber. 


The compressive stress in the composite unit cell in Y direction due to q can be expressed 
by (J 2 = j= E° n £ 2 = E° 22 - E\ 2 ~ (2-39) 


where a 2 is transverse stress due to q, £ is transverse strain due to q, T is the thickness 
of the composite, and E 22 is transverse modulus of the composite. 


Substituting q from (2-39) to (2-38), we obtain 

8 C E°T8C 

-r<j f A f =- — 


(2-40) 


and 
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8 C=- 


(2-41) 


8C 2 a f A f 

e° 22 i 2 t 


Assuming the composite and the fiber bundle have the same strain in the X direction. 



If fiber volume fraction is denoted as Vf , then 

A f =A 0 V f =CTV f (2-42) 

Therefore, (2-34) can be rewritten as 


8C=- 


8C 3 E f (7 0 V f 

p° ] 2 


(2-43) 


Effective strain £, can be expressed as 


£. = 


izl o 
lo 
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^l+F(l,C)8C 
V E\\ j 

f a , 8 C 3 E f o^V f F(l,C) 
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— n 1 1 — r\ — n 




) 


(2-44) 


Therefore 



(2-45) 


where AT is a knockdown factor of longitudinal modulus for the composite. 
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(2-46) 


8 C 3 E f V f F(l,C) 

1 + — - 

E°l 3 
22 l 

Equation (2-46) can be used to estimate the waviness knockdown factor for misoriented 
unidirectional composite segments with axially loaded wavy tows. The value of K is less 
than one and depends on the ratios of C/l and E f j E° 2 . The larger the C// and/or 

E f / E , the smaller the K. 
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Chapter 3 


Prediction of 3-D 
Woven and Braided 
Composite Properties 

This chapter outlines a preliminary implementation of the 3-D composite mechanics 
module. Elastic properties and stress limits of three-dimensionally reinforced fiber 
composites are evaluated. Fiber orientations of weave and braid patterns are defined with 
reference to composite structural coordinates. Orthotropic ply properties and stress limits 
computed via micromechanics are transformed to composite structural coordinates and 
integrated to obtain the 3-D properties. Various composites and fiber reinforcement 
patterns are considered. Results are compared with experimental data from the literature. 

3.1 Methodology 

A composite mechanics module was implemented to compute the elastic properties and 
stress limits of general 3D fiber composites. The implementation uses micromechanics 
modules to evaluate inlraply hybrid composite properties. The intraply properties are 
computed separately for each fiber orientation present in the 3D braided/woven 
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composite. The assembly of 3D braided composite properties is accomplished by rotating 
for each braid the stress-strain relations from the local coordinate axis of the braid into the 
structural coordinate system of the 3D composite via tensor transformations. The local 
coordinate axis of a braid is taken with the x axis along the fiber direction of the braid. 
For each braid orientation, the computational module input file is modified to define the 
direction of the braid x axis by inputting the angles made by the braid axis with the three 
global composite structural coordinates. Unless the local x axis is in the global z direction, 
the local y axis of the braid is determined in the direction of the cross product of the global 
z axis and the local x axis. If the local x-axis is in the global z direction, the local y axis is 
taken to be in the direction of the cross product of the local x axis by the global x axis. 
The local z axis is determined by taking the cross product of local x and local y directions. 
Components of the unit vectors of braid local coordinates written with respect to the 
composite structural coordinates define the coefficients of the braid orientation matrix, see 
chapter 2 for details. 

The properties of each braid are transformed to the structural directions as described in 
the previous chapter, multiplied by the ratio of fibers in that braid to the total composite 
fibers, and superimposed onto the 3D composite properties. When the contribution of 
each braid is added, composite structural properties are obtained. Fiber volume ratios may 
be different for each braid orientation as specified by different material cards in the input 
file. Also, different fibers may be selected for different braids and laminae. Results 
indicate the implemented method is consistent with published results [Cox and Dadkhah, 
1995]. There remain some differences between the computed elastic properties and 
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experimental data. These may be attributed to the effect of residual stresses from the 


curing process as well as the waviness of fibers. 


The fundamental distinction of the present method is that the augmented composite 
mechanics code is also able to compute the 3D composite stress limits by accumulating the 
contribution of each braid/weave/ply to each composite strength component via tensor 
transformations of local strengths in the absolute value. The stress limits in the composite 
structural coordinates corresponding to first ply failure (FPF) are written in general as: 


7=1 


Where Sc± are the positive or negative stress limits in the composite structural coordinate 
directions; S w± are the stress limits in the local coordinate system of each braid or ply; and 



is the property rotation matrix taken in absolute value. In composite strength 


computations, the fiber and matrix contributions are considered separately to enable the 
representation of matrix damage that degrades the composite properties prior to fracture. 
Each of the three composite normal stress limits are computed in both tension and 
compression. At the stage of matrix damage initiation, tensile stress limits consider the 
contribution of all components except for fiber tensile strength. On the other hand, in 
computing the composite tensile fracture strengths, only the contribution of fiber strengths 
are considered. Compressive normal stress limits are also computed both with reference 
to matrix crushing as well as with regard to fiber fracture. Similarly, shear strengths with 
reference to the three global structural axes are obtained by the accumulation of absolute 
tensor transformations of braid shear strengths as well as the relevant normal strength 
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contributions. In the computation of shear strengths, the appropriate normal strength 
components are the compressive strength in the braid longitudinal directions and tensile 
strengths in the braid transverse directions. For each one of the eleven specimens 
considered, the computed elastic properties and the composite stress limits are tabulated in 
Appendix B. The computed components of strength the names of which end with the 
letter "f correspond to the stress limits for fiber fracture at the absence of matrix 
degradation in the composite. On the other hand, strength components that are assigned 
names ending with the letter "m" indicate the matrix damage initiation stresses. Any 
general 3D composite structure with arbitrary fiber orientations can be analyzed for stress 
limits as well as composite elastic properties. For the woven composite examples 
considered in this chapter, the x axis is taken in the direction of the stuffer fibers, the y axis 
is taken in the direction of the filler fibers, and the z axis is taken in the normal direction. 
Warp fibers are in various different directions depending on the specimen considered. 

3.2 Results and Discussion 

The developed method was tested by simulating an experimental study on eleven woven 
composite structures (Cox and Dadkhah 1995). Of the eleven specimens five were lightly 
compacted with fiber volume ratio of approximately 0.35 whereas the remaining six were 
heavily compacted with fiber volume ratios of approximately 0.60. The specimens were 
designated as indicated and described in Table 3.1. 
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Table 3. 1 Designation of Woven Composite Types 


SYMBOL DESCRIPTION 


LL 

LT 

LO 

HL 

HT 

HO 


Lightly Compacted Layer-to-layer Angle Interlock 
Lightly Compacted Through-the-thickness Angle Interlock 
Lightly Compacted Through-the-thickness Orthogonal Interlock 
Heavily Compacted Layer-to-layer Angle Interlock 
Heavily Compacted Through-the-thickness Angle Interlock 
Heavily Compacted Through-the-thickness Orthogonal Interlock 


Table 3.2 shows the fiber volume ratios for the different tows of each specimen. Appendix 
A shows the properties of fibers and matrices which were used in this simulation. The 
stuffer and filler fibers were made of AS graphite for all specimens. The LL and LT 
woven composite types each contained two specimens of which one had warp weaver 
fibers made of graphite and the other had warp weaver fibers made of glass. All fibers 
were made of graphite for the heavily compacted composites. Each type of heavily 
compacted weave pattern had two specimens with different ratios of the stuffer/filler/warp 
weaver fiber ratios. 

The implemented composite mechanics module was able to predict 3D composite 
properties and stress limits consistently. The computed results for all eleven specimens are 
presented in tabular form in Appendix B. The present results are compared with those 
computed by Cox and Dadkhah and those obtained from experiments. In the majority of 
cases the new module evaluated the 3D stiffness properties closer to experimental 
measurements compared to the orientation averaging method reported by Cox and 
Dadkhah. However, in the direction of filler yarns, stiffness computations by the extended 


NASA/CR— 2005-213963 


43 



code, as well as those computed by Cox and Dadkhah (1995), were considerably higher 
than the experimentally measured values. Composite stress limits based on fiber strength 
and matrix strength computed by the composite mechanics code were always able to 
provide upper and lower bounds, respectively, to the experimental failure strengths. 

The availability of the present 3D composite analysis method enables the assessment of 
damage tolerance, as well as structural response evaluation of braided and woven 
composites. However, the developed methods can be further improved by quantifying the 
effects of tow waviness on composite structural response and stress limits. Also, the 
alternative formulation proposed in the previous chapter that preserves the spatial 
configurations of individual stuffer and filler tows may be more appropriate for composites 
subjected to bending. 
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Table 3.2 Fiber Volume Fraction of Specimens Considered 



Specimen 

Thickness 

Nominal Volume 
Fractions 

Total Fiber 
Volume 
Fraction e 

Fraction by Volume of all 
Fibers Lying in Stuffers , 
Fillers, Wrap Weavers f 

Composite 

(cm) a 

gm 

v f c 

M 

V 

fs 

fr 

fw 

1-L-l 

1.26 

0.14 

0.15 

0.07 

0.35±0.03 

0.385 

0.418 

0.197 

l-L-2 

1.24 

0.14 

0.20 

0.05 

0.370±0.005 s 

0.066±0.004 h 

0.347 

0.502 

0.151 

1-T-l 

1.02 

0.16 

0.21 

0.05 

0.466±0.003 

0.381 

0.504 

0.115 

l-T-2 

0.97 

■ 

0.22 

0.04 

0.408±0.020 g 

0.044±0.004 h 

0.406 

0.497 

0.097 

1-0 

0.88 


0.25 

0.04 

0.483+0.010 

0.387 

0.524 

0.090 

h-L-1 

0.561 

0.38 

0.22 

0.05 

0.620+0.008 

0.587 

0.340 

0.073 

h-L-2 

0.625 

0.33 

0.21 

■ 

0.557+0.015 

0.580 

0.375 

0.045 

h-T-1 

0.573 

0.37 

0.22 

0.065 

0.613+0.003 

0.571 

0.331 

0.098 

h-T-2 

0.577 

0.36 

0.23 

0.035 

0.592+0.014 

0.571 

0.369 

0.059 

h-O-1 

0.579 

0.37 

0.22 

0.045 

0.619+0.008 

0.586 

0.340 

0.073 

h-O-2 

0.587 

0.35 

0.23 

0.065 

0.593+0.014 

0.545 

0.353 

0.102 


• a in direction normal to wrap and weft direction 

• b V s = volume of fraction of stuffer (straight wrap) twos 

• c V f = volume of fraction of filler (weft) tows 

• d V w = volume fraction of wrap weaver (3D warp) tows 

• e measured by acid digestion 

• f determined from weaver's specification 

• g graphite fibers 

• h glass fibers 


NASA/CR— 2005-213963 


45 



































































































Chapter 4 


Damage Propagation in 3-D 
Woven Composite Structures 


The design evaluation of composite structures requires an assessment of their safety and 
durability under service loads and possible overload conditions. The ability of designing 
composites with numerous possible fiber orientation patterns, choices of constituent 
material combinations, tow/ply drops and hybridizations, render a large number of possible 
design parameters that may be varied for an optimal design. The multiplicity of composite 
design options presents a logistical problem, prolonging the design and certification 
process and adding to the cost of the final product. It is difficult to evaluate composite 
structures due to the complexities in predicting their overall congruity and performance, 
especially when structural degradation and damage propagation take place. The 
predictions of damage initiation, damage growth, and propagation of fracture are 
important in evaluating the load carrying capacity, damage tolerance, safety, and reliability 
of composite structures. Quantification of the structural fracture resistance is also 
fundamental for evaluating the durability/life of composite structures. The most effective 
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way to obtain this quantification is through integrated computer codes that couple 
composite mechanics with structural analysis and damage progression modeling. 

An important feature of computational simulation is the assessment of damage stability or 
damage tolerance of a structure under loading. At any stage of damage progression, if 
there is a high level of structural resistance to damage progression under the service 
loading, the structure is stable with regard to fracture. The corresponding state of 
structural damage is referred to as stable damage. On the other hand, if damage 
progression does not encounter significant structural resistance, it corresponds to an 
unstable damage state. Unstable damage progression is characterized by very large 
increases in the amount of damage due to small increases in loading. Whereas during 
stable damage progression the amount of increase in damage is consistent with the 
increase in loading. 

Internal damage in composites is often initiated as cracking due to normal stresses 
transverse to fiber orientation. Further degradation is in the form of additional cracking, 
delaminations, and fiber fractures that may lead to structural fracture. Because of the 
numerous possibilities with material combinations, composite geometry, fiber orientations, 
and loading conditions, it is essential to have an effective computational capability to 
predict the behavior of composite structures for any loading, geometry, composite 
material combinations, and boundary conditions. The predictions of damage initiation, 
growth, accumulation, and propagation to fracture are important in evaluating the load 
carrying capacity and reliability of composite structures. Quantification of the structural 
fracture resistance is also required to evaluate the durability/life of composite structures. 
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The behavior of laminated composite structures under various loading condition is studied 
in Minnetyan et al. (1990-1995). The proposed research will extend the scope of a 
computational tool that has been developed to examine the progressive damage response 
of conventional laminated fiber composites to the evaluation of 3-D woven composite 
structures. This integrated computational tool can effectively be broken into three main 
modules: 

1. A micromechanics and macromechanics module of 3D composites. 

2. A structural analysis module. 

3. A damage progression tracking module. 

Details of the first module are discussed in the previous chapter. This chapter discusses the 
second and third modules. 


4.1 Structural Analysis Module 

The main tool of the structural analysis module is the finite element method. After the 
micromechanics analysis module generates the elastic properties for 3-D woven 
composite, a finite element analysis module is called to analyze the entire structural 
response based on the composite properties. In general, the type of finite element model 
used depends on the complexity of the structure and the availability of computer 
resources. There are two possible choices for the analysis of composite structures. One is 
using anisotropic 3-D solid elements (such as hexahedral or brick elements) that accept the 
computed 3-D composite properties directly. However, the modeling of a practical 
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composite structure with three dimensional finite elements requires huge computer 
memory and CPU time and it is usually impractical for analyzing real structures. The 
second option is to use anisotropic shell elements that use the composite plate/shell 
element properties obtained in the previous chapter. Anisotropic plate or shell elements 
represent through-the-thickness properties of the three-dimensional composite and 
therefore are much more efficient compared to three dimensional elements. Therefore, 
implementation is focused on the use of plate/shell elements in the finite element analysis 
module. The finite element module accepts the 3-D woven composite force-deformation 
relations predicted by the composite macromechanics module, and carries out a stress 
analysis to generate the generalized stresses Nx, Ny, Nxy, Mx, My, Mxy for each node. 
The generalized stresses are supplied back to the composite mechanics module for the 
computation of local stresses in the primary and weaver domains and failure analysis. 

4.2 Failure Criteria 

Progressive damage and fracture simulations will be carried out by imposing failure 
criteria locally within each micromechanics subvolume. Micromechanics subvolumes will 
be obtained by subdividing each micromechanics volume into regions with characteristic 
fiber configuration. Within each subvolume local coordinate orientation in the material 
directions will be identified. At each individual load step, the stuffer, filler, and weaver 
stresses and strains obtained through 3D woven composite microstress analysis are to be 
checked according to distinct failure criteria. The first twelve failure modes are associated 
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with the positive and negative limits of the six local stress components in the material 
direction as follows: 

■^nic < °)n < S/ut 
S(22 C < ' ®122 ^ $122T 
S( 33 c < < S l33T 

^m(-) < < ^i2(+> 

^ 23 (-) ^ °/23 < S(2X+) 

S(l3(-) < *3(13 < ^ 13 (+) 

The thirteenth failure mode is a combined stress failure criterion, or modified distortion 
energy (MDE) failure criterion which is obtained by modifying the usual distortion energy 
failure criterion. The modification takes into account the significant differences in the 
stress limits of the longitudinal and transverse directions of an orthotropic composite ply. 
Each component of ply stress is normalized with respect to its limiting strength. No 
relationship is assumed between normal and shear strengths. The MDE failure criterion 
has been demonstrated to be a good predictor of combined stress failure in composites. It 
may be considered as a variation of the Tsai-Hill theory (Tsai 1968, Hill 1950). The MDE 
failure criterion (Chamis 1969) can be expressed as: 
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where a and (5 indicate tensile or compressive stresses, S (lloc is the local longitudinal 
strength in tension or compression, S (22a is the transverse strength in tension or 
compression, and the directional interaction factor is defined as: 

— v' (l + 4 u 12 — V\i)E 22 + (l ~ v 2 ^)E u 

~ [E»E a (2+v a + 1>„)(2 + 1> 2 , + uj ] 1/2 

where K' lna p is a theory-experiment correlation factor. 

The directional interaction factor reduces to unity for homogeneous isotropic materials. 

4.3 Simulation of Damage Progression 

After each finite element stress analysis, failure criteria will be used to evaluate possible 
failure at each primary and weaver domain of each ply at each node of the composite 
structure. Once the damage at each node is assessed, a damage index will be created to 
record the damage information for each damaged node. The damage index will contain the 
node number, the ply number, and the list of damage criteria that have become activated. 
When a new failure occurs at certain stuffer, filler, or weaver domain, the damage index 
will be updated correspondingly. The properties of each domain will be degraded 
according to their damage index. 

If there is no damage after a load increment, the structure will be considered to be in 
equilibrium and an additional load increment will be applied. Figure 4. 1 shows a schematic 
of damage tracking, expressed in terms of a load-displacement relationship. Point 1 
represents the last equilibrium state before initial damage. When the structure is loaded by 
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an additional load increment to point 2, ply failure criteria indicate damage initiation. The 
composite properties affected by the damage are degraded, the computational model is 
reconstituted with updated finite element mesh and material properties, and the structure is 
reanalyzed under the same load increment to reach point 3. However, at point 3, 
composite ply criteria indicate additional damage. Accordingly, structural properties are 
further degraded and analysis is repeated under the same load increment to reach point 4. 
There is no further damage at point 4 and the structure is now considered to be in 
equilibrium with the external loads. Since the structure is now in equilibrium, a new load 
increment is applied to reach point 5. At this point, the process of checking the failure 
criteria, updating the finite element mesh, and updating the material properties is repeated 
with possible damage growth and propagation. 



Figure 4. 1 Schematic of damage tracking of woven composites 
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The damage progression module keeps a detailed account of composite degradation for 
the entire structure and also acts as the master executive module that directs the 
composite mechanics module to perform micromechanics, macromechanics, laminate 
analysis and synthesis functions. The damage progression module also calls the finite 
element analysis module with thick shell analysis capability to model laminated composites 
for global structural response. A schematic of the simulation cycle is shown in Figure 4.2. 

4.4 Damage Energy Release Rate 

The measure of global fracture toughness is defined as a Damage Energy Release Rate 
(DERR) that is equal to the incremental amount of damage energy expanded for the 
creation of unit damage volume in the composite structure. The magnitude of the DERR 
varies during progressive degradation of the composite structure under loading, reflecting 
the changes in the fracture toughness of the composite. Computation of DERR during 
progressive fracture is needed in order to evaluate the composite fracture toughness and 
the degree of imminence of failure. 

The DERR is computed as the ratio of incremental damage energy to the corresponding 
incremental damage volume that is generated. 

DERR - Dama8eEner8y 
Damage Volume 

Damage energy can be computed as follows: 

Damage Energy = £ X (° 5S) / Ej )v, 
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Where Sj is a local composite strength, E s is modulus corresponding to S p Vi is the volume 
of damage. The traverse of a local minimum value by the DERR during the progression of 
fracture typically precedes a very high rate of damage propagation and generally predicts 
the imminence of total failure. Typically, at the stage of damage initiation, there is a high 
rate of energy release that dissipates a significant portion of the strain energy stored in the 
composite structure before the onset of initial damage. After the first burst of energy 
release, DERR usually drops down to a lower level, indicating the significant reduction in 
fracture stability of the damaged composite under the applied loading. The DERR usually 
reaches its peak value as global failure occurs. 
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Check if the entire structure 
is collapsed 


Print out the ultimate 
loads and all relative 
information 


Separate all the element 
connecting to the fractured 
nodes, created additional 
nodes and renumber the nodes 
and element connectivities 


Degrade the material properties 
for the damaged nodes according 
to the damage model 


Figure 4.2 Flow chart of integrated computer code for simulating 
damage propagation of 3D woven composites 
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Chapter 5 


DAMAGE PROGRESSION 
IN CARBON-FIBER 
PLASTIC I-BEAMS 

5.1 Introduction 

For a rational design, it is necessary to quantify the defect and damage tolerance of a 
candidate structure. The assessment of damage tolerance requires a capability to simulate 
the progressive damage and fracture characteristics of structures under loading. Critical 
components of a structure are required to remain safe and to be able to function under 
loading after experiencing some damage. The cause of damage may be an accident, defect, 
or unexpected overloading. Damage tolerance of a structure is quantified by the residual 
strength; that is, the additional load carrying ability after damage. Composite structures 
are well suited for design with emphasis on damage tolerance because continuous fiber 
composites have the ability to arrest cracks and prevent self-similar crack propagation. For 
most fiber reinforcement configurations, cracks and other stress concentrators do not have 
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as important an influence on composites as they do for homogeneous materials. Another 
important factor is the multiplicity of design options for composites. The ability of 
designing composites with numerous possible fiber orientation patterns, choices of 
constituent material combinations, ply drops, and hybridizations render a large number of 
possible design parameters that may be varied for an optimal design. 

Graphite/epoxy composites are used in various retrofit and original construction projects, 
as well as other structural applications such as aircraft structures, automotive vehicles, 
pressurized containers, special purpose enclosure structures, and so on. In most 
applications it is important to achieve low weight, high strength, adequate stiffness, and 
safety when composite structures are required to withstand significant loads. In the current 
study, discussion is focused on composite beams formed into a structural shape subject to 
three-point bending. Damage initiation, growth, accumulation, propagation to fracture, 
and buckling stability are evaluated. The damage initiation load and the maximum 
structural fracture load or the buckling instability load are quantified for three different 
laminate configurations. 

Design considerations with regard to the progressive fracture of composite structures 
require an evaluation of damage initiation and propagation mechanisms under expected 
service loading. Concerns for safety and survivability of critical components require a 
quantification of the structural fracture resistance under sustained loads. A significant 
design parameter with regard to composite damage tolerance is the laminate configuration. 
In general, quasi-isotropic laminates yield better damage tolerance. However, in many 
cases, a quasi-isotropic laminate may not be the most efficient with regard to structural 
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strength and performance when there is no damage. For a rational design process, it is 
necessary to quantify the structural damage tolerance for a candidate design. Damage 
initiation and progression characteristics are much more complex for laminated fiber 
composites compared to homogeneous or orthotropic materials. The structural fracture 
process of a laminated composite depends on many parameters such as laminate 
configuration, fiber volume ratio, constituent stifftiess/strength/hygrothermal parameters, 
and the fabrication process. Recent developments in computational simulation technology 
have made it possible to evaluate the details of progressive damage and fracture in 
composite structures. Computational simulation enables assessment of the damage 
initiation and propagation loads. A damage energy release rate is evaluated globally during 
simulation by computing the work done by the applied loads per unit damage created. The 
damage energy release rate is used to quantify the structural damage tolerance at different 
stages of degradation. The influence of local defects or flaws, through-the-thickness 
cracks, and effects of the fabrication process in terms of residual stresses can be taken into 
account. 

Laminated composite design practice has been based on extensive testing with attempts to 
apply formal fracture mechanics concepts to interpret test results. In certain cases, 
interpretation of composite test data via fracture mechanics has been satisfactory. 
However, in most cases, fracture mechanics methods have significantly mispredicted the 
strength of fiber composites. Reconciliation of test results with fracture mechanics has 
required significant modifications of effective fracture toughness and the calibrations of 
specific, laminate configuration dependent, effective stress concentration field parameters. 
Additionally, required adjustments of fracture mechanics parameters have had to be 
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reassessed with every change in constituent and laminate characteristics. The complete 
evaluation of laminated composite fracture requires an assessment of ply and subply level 
damage/fracture processes. 

The present approach bypasses traditional fracture mechanics to provide an alternative 
evaluation method, conveying to the design engineer a detailed description of damage 
initiation, growth, accumulation, and propagation that would take place in the process of 
ultimate fracture of a fiber composite structure. Results show in detail the damage 
progression sequence and structural fracture resistance during different degradation 
stages. This study demonstrates that computational simulation, with the use of established 
material modeling and finite element modules, adequately tracks the damage growth and 
subsequent propagation to fracture for carbon-fiber reinforced plastic (CFRP) composite 
I-beams subjected to three point bending. 


5.2 Methodology 

An integrated computer code has been developed to simulate damage initiation, damage 
growth, and fracture in composite structures under various loading and environmental 
conditions (Chamis et al. 1996). Computational simulation has been used for the 
investigation of the effects of composite degradation on structural response (Minnetyan et 
al. 1991), composite structures' global fracture toughness (Minnetyan et al. 1990), the 
effect of the hygrothermal environment on durability (Minnetyan et al. 1992b), damage 
progression in composite shells subjected to internal pressure (Minnetyan et al. 1992a), 
damage of composite cylinders subjected to mechanical loads as well as internal pressure 
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(Gotsis et al. 1996), the behavior of discontinuously stiffened composite panels under 
compressive loading (Minnetyan et al. 1995), and damage propagation in thick composite 
shells under external pressure (Minnetyan and Chamis, 1997). The present paper examines 


damage and fracture progression in CFRP I-beams with three different layups. 



Figure 5.1 Computational Simulation Cycle 

The developed computer method consists of an integrated code with three modules: (1) 
composite mechanics; (2) finite element analysis; and (3) damage progression modeling. 
The overall evaluation of composite structural durability is carried out in the damage 
progression module that keeps track of composite degradation for the entire structure. 
The damage progression module relies on a composite mechanics code (Murthy and 
Chamis, 1986) for composite micromechanics, macromechanics, and laminate analysis, 
and calls a finite element analysis module that uses anisotropic thick shell elements to 


NASA/CR— 2005-213963 


60 











model laminated composites. 


Figure 5.1 shows a schematic of the computational simulation cycle. The composite 
mechanics module is called before and after each finite element analysis. Prior to each 
finite element analysis, the composite mechanics module computes the composite 
properties from the fiber and matrix constituent characteristics and the composite layup. 
The finite element analysis module accepts the composite properties that are computed by 
the composite mechanics code at each node and performs the analysis at each load 
increment. After an incremental finite element analysis, the computed generalized nodal 
force resultants and deformations are supplied to the composite analysis module that 
evaluates the nature and amount of local damage, if any, in the plies of the composite 
laminate. Individual ply failure modes are assessed by using failure criteria associated with 
the negative and positive limits of the six ply-stress components in the material directions. 
In addition to the failure criteria based on stress limits, interply delamination due to 
relative rotation of the plies and a modified distortion energy (MDE) failure criterion that 
takes into account combined stresses are considered (Murthy and Chamis, 1986). 
Depending on the dominant term in the MDE failure criterion, fiber failure and/or matrix 
failure is assigned. 

If the dominant term in the MDE failure criterion is ply longitudinal tensile or compressive 
stress, both fiber and matrix failures are simulated, and the ply stiffness is lowered to a 
negligible value. On the other hand, if the dominant term in the MDE criterion does not 
correspond to a longitudinal tensile or compressive failure mode, only the matrix stiffness 
is reduced for the damaged ply, and the fiber stiffness is retained. A ply that sustains 
matrix damage is able to carry longitudinal tensile stresses but would fracture 
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Loading 



Figure 5.2a CFRP I-Beam and Loading 


outside flange 



Figure 5.2b Cross Section of CFRP I-Beam 
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Duplicate nodes 



when subjected to longitudinal compression. The generalized stress-strain relationships are 
revised locally according to the composite damage evaluated after each finite element 
analysis. The model is automatically updated with a new finite element mesh having 
reconstituted properties, and the structure is reanalyzed for further deformation and 
damage. If there is no damage after a load increment, the structure is considered to be in 
equilibrium, and an additional load increment is applied, leading to possible damage 
growth, accumulation, or propagation. Simulation is continued until global structural 
fracture. 

During progressive damage tracking, the following terminology is utilized to describe the 
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various stages of degradation in the composite structure: (1) damage initiation refers to 
the start of damage induced by loading that the composite structure is designed to carry; 
(2) damage growth is the progression of damage from the location of damage initiation to 
other regions; (3) damage accumulation is the increase in the amount of damage in the 
damaged regions with additional damage modes becoming active; (4) damage propagation 
is the rapid progression of damage to other regions of the structure; and (5) structural 
fracture is the ultimate disintegration of the specimen. 

5.3 Carbon-Fiber Reinforced Plastic I-Beams 

Computational simulation examples in the present paper consist of composite I-beams for 
which experimental results are available from the published literature. Three types of I- 
beams from Takayanagi et al. (1994) were simulated. The reported I-beams were 
fabricated from Toray F3060 (T300/ 3601) prepreg cured at 180°C. The stacking 
sequences of the I-beams are shown in Table 5.1, and a schematic view of an I-beam is 
shown in Figure 5.2. The layup for the web and the inside flange was a ±45° angle-ply for 
all three specimens. The number of plies for the web was 16. The number of plies for the 
inside flange was eight in all cases. Three different layup schemes for the outside flange 
were considered. The layup for specimen A was unidirectional (UD), the layup for 
specimen B was quasi-isotropic (QI), and the layup for specimen C was angle plied (AP) 
(see Table 5.1 for the precise definitions of UD, QI, and AP laminate configurations). 
Table 5.2 shows the dimensions for the cross sections of the I-beams with reference to 
Figure 5.2(b). 


NASA/CR— 2005-213963 


64 



Table 5.1 Laminate Stacking Sequence of the CFRP I-Beams 


Specimen 

Flange 

Web 

Outside 

Inside 

A 

UD: [0] 8 

AP:[+45/-45] 4 

AP:[(±45) 4 ] s 

B 

QI: [45/0/-45/90] 2 

AP:[+45/-45] 4 

AP:[(±45) 4 ] s 

C 

AP: [+45/-45] 4 

AP:[+45/-45] 4 

AP:[(±45) 4 ] s 


Table 5.2 Dimensions of the Cross-Sections of the CFRP I-Beams 


Specimen 

b (mm) 

b w (mm) 

t f (mm) 

t w (mm) 

A 

27.56 

2.48 

2.10 

60.17 

B 

28.32 

2.24 

2.00 

60.17 

C 

28.36 

2.46 

2.06 

60.17 


A finite element model of the I-beam is shown in Figure 5.3. Mindlin type thick composite 
shell elements were utilized to represent the through- the-thickness laminates. The upper 
and lower flanges were uniformly divided into 200 elements each. The web was divided 
into 250 uniform sized elements. Duplicate nodes were used to represent the connection 
of the web to the flanges. Figure 5.3 shows the duplicate nodes at a typical section of the 
I-beam. In the computational model, each pair of duplicate nodes is described by assigning 
a slave-master relationship that is implemented by eliminating the degrees of freedom of 
the slave node from the stiffness equations and adding the stiffness contributions from the 
slave node to those of the master node. During computational simulation of damage 
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progression, if one of the nodes of a duplicate node pair was fractured, the associated 
duplicate node relationship was terminated. 

5.4 Fiber and Matrix Properties Calibration 

To enable constituent level damage tracking in the computational simulation, it was 
necessary to identify the matrix/fiber constituent properties and stress limits. Takayanagi et 
al. (1994) gave the ply properties, but did not provide explicit information on fiber and 
matrix properties. The composite mechanics module of this code (Murthy and Chamis 
1986) was used to calibrate the fiber and matrix properties so that the computed ply 
properties matched those of the T300/3601 composite given by Takayanagi et al. (1994). 
The resident databank of the composite mechanics code provided properties of similar 
fibers and matrices. Properties of T300 from the databank were used as the initial values 
of fiber T300, which was used in this study; properties of matrix 3501 (since matrix 3601 
was not included in the ICAN databank) were used as initial values for matrix 3601. Using 
the properties from the databank as initial trial values, ply properties and composite 
properties were identified via a least squares error sum minimization method. 

Ply properties computed by the composite mechanics module were compared with those 
experimentally observed by Takayanagi et al. (1994). If the differences between the 
computed and observed ply properties was greater than acceptable, the fiber and matrix 
properties were revised, and ply properties were recomputed until the actual ply properties 
were matched with sufficient accuracy. However, because of the natural complexity of the 
composite material, as well as the uncertainty and variability of experimental observations. 
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it was not possible to obtain computed ply properties that were exactly the same as those 
experimentally observed. In the present study, the difference between the computed 
properties and experimental observations was minimized in the root-mean-square relative 
error sense. The root-mean-square relative error Q was defined as 




rj observed ^ computed 

i * i 


observed 

n ,=i 

k 

* J 


where p i observed = value of composite property (such as En, E 22 , Gi 2 , etc.) from experiment; 
p computed _ va j ue Q f same property computed by composite mechanics module; and n = 
number of composite ply properties observed. 

The smaller the value of Q, the better the agreement between the computed and observed 
ply properties, and therefore more accurately identified are the fiber and matrix constituent 
properties. The fiber and matrix properties identified by minimizing Q are given in Tables 
5.3 and 5.4, respectively. 

Based on the fiber and matrix properties given in Tables 5.3 and 5.4, the composite 
properties were computed by the composite mechanics code. Tables 5.5, 5.6, and 5.7 
show comparisons of computed and observed composite properties as well as the 
associated percent relative error Q. 
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Table 5.3: Fiber Properties 


Parameter 


Value 


Number of fibers per end 
Fiber diameter 
Fiber Density 

Longitudinal normal modulus 
Transverse normal modulus 
Poisson's ratio Vn 
Poisson's ratio V 23 
Shear modulus Gn 
Shear modulus G 23 

Longitudinal thermal expansion coefficient 
Transverse thermal expansion coefficient 
Longitudinal heat conductivity 
Transverse heat conductivity 
Heat capacity 
Tensile capacity 

Compressive strength 


3000 

0.00762 mm(0.300E-3 in) 

1771 kg/m (0.064lb/in ) 

197 GPa(28.5E+6 psi) 

17.2 GPa(2.5E+6 psi) 

0.20 

0.25 

8.96 Gpa ( 1.3E+6 Psi ) 

4.82 Gpa ( 0.7E+6 Psi ) 

-0.99E-6°C ( -0.55E-6°F ) 

1.01E-6°C( 0.56E-5°F ) 

43.4 J-m/hr/m /°C (580BTU-in/hr/in/°F ) 
4.34 J-m/hr/m 2 /°C ( 58BTU-in/hr/in/°F ) 
712 J/kg/°C (0.17BTU/ lb/°F ) 

2.41 Gpa ( 350 ksi) 

2.07 GPa ( 300 ksi ) 


able 5.4: Matrix properties: 


Parameter 

Value 

Matrix density 

1265 kg/m 3 (0.0457 Ib/in) 

Normal modulus 

3.79 GPa ~(550ksi) 

Poisson's ratio 

0.34 

Coefficient of thermal expansion 

0.72E-4/°C (0.4E-4/°F ) 

Heat conductivity 

1 .2 5 BTU- in/hr/in 2 /°F 

Heat capacity 

0.25 BTU/lb/°F 

Tensile strength 

114 MPa (16.5 ksi) 

Compressive strength 

423 MPa(61.3ksi) 

Shear strength 

148 MPa (21.4ksi) 

Allowable tensile strain 

0.02 

Allowable compressive strain 

0.05 

Allowable shear strain 

0.04 

Allowable torsional 

0.04 

Void conductivity 

16.8 J-m/hr/m 2 / °C (0.225BTU-in/hr/in/°F) 

Glass transition temperature 

216°C (420°F ) 
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Table 5.5 Elastic Constants of the Laminates 


Type of Laminate 

Eu(GPa) 

E 22 (GPa ) 

G I2 (GPa) 

Q(%) 

V f (%) 

UUD 

Experimental 

100.1 

9.12 

3.71 

5.6 

49.8 


Simulated 

99.77 

8.43 

3.49 



QQI 

Experimental 

40.6 

40.6 

14.9 

4.1 

49.4 


Simulated 

38.6 

38.6 

14.7 



AAD 

Experimental 

11.8 

11.8 

26.1 

4.8 

50.3 


Simulated 

12.5 

12.5 

26.3 




Table 5.6 Elastic Constants of the Flange Laminates 


Specimen 

E x (GPa) 

E y (GPa) 

G„(GPa) 

Q(%) 

V f (%) 

AA 

Experimental 

72.3 

13.1 

18.6 

2.6 

62.8 


Simulated 

71.2 

12.6 

18.6 



BB 

Experimental 

35.5 

34.8 

27.0 

2.6 

66.0 



35.9 

35.9 

26.8 



CC 

Experimental 

15.5 

15.5 

34.3 

2.6 

66.1 


Simulated 

16.0 

16.0 

34.3 




Table 5.7 Elastic Constants of the Web Laminates 


Specimen 

E w (GPa) 

G^GPa) 

Q(%) 

V f (%) 


Experimental 

11.7 

25.9 

4.24 

50.0 

A 

Simulated 

12.4 

26.1 




Experimental 

13.8 

30.5 

2.0 

58.7 

B 

Simulated 

14.2 

30.5 




Experimental 

12.3 

27.1 

3.5 

52.3 

C 

Simulated 

12.9 

27.3 
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5.5 Results and Analysis 


Damage initiation and progression were monitored as the composite I-beam specimens 
were loaded at midspan. The loading configuration was as shown in Figure 5.2(a). In the 
computer model, the knife-edge supports were described as follows. For the left end 
support, bottom flange nodes in contact with the knife-edge support were restrained 
against translation in the plane of the web (pin supports). Additionally, the bottom flange 
node immediately under the web was also restrained against translation perpendicular to 
the plane of the web. For the right end support, all bottom flange nodes in contact with the 
knife-edge were restrained in the load direction (roller supports). The loading nose was 
represented by an attached finite element model of a high-strength steel plate that was 20 
mm (0.787 in.) thick and 41 mm (1.6 in.) high. The attachment of the steel plate finite 
elements to the top flange of the beam was described using duplicate node constraints. 
The load was applied to the top of the loading nose. It was started at 0.20 kN (46 lb) and 
increased gradually. Specimen A was simulated first. When the load applied on specimen 
A reached 2.56 kN (575 lb), damage initiation occurred at the midspan of the beam, in the 
upper flange by matrix tensile failures due to the transverse tensile fracture of the 0° 
outside plies immediately above the web. Damage grew to adjacent nodes in the transverse 
direction as the load was increased to 2.74 kN (616 lb). However, the damage mode 
remained with matrix cracking in the 0° outside plies only. When the load was increased to 
6.36 kN (1,430 lb), the upper flange outermost ply under the loading nose developed fiber 
fractures in the longitudinal compressive failure mode. At 6.81 kN (1,530 lb), ply 
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compressive fiber fractures progressed to the three outermost 0° plies of the top flange. 
Also, fiber tensile fractures began to occur at the bottom flange at midspan, immediately 
under the web, in the outermost 0° plies. At 7.65 kN (1,720 lb), all eight plies of the outer 
top flange experienced compressive failures under the load. At the bottom flange, tensile 
failures grew to the outermost four plies. At 7.96 kN (1,790 lb) loading, ply compressive 
fractures grew into through-the-thickness laminate cracks of the upper flange at midspan. 
First, laminate cracking began at the edges of the top flange; second, the center flange 
under the loading nose developed cracking; and third, the intermediate nodes between the 
center and the edges of the top flange were cracked. In the bottom flange, fiber tensile 
failures progressed across the flange width in the outermost three plies. Damage in the 
center nodes of the bottom flange remained the same with four fractured plies. As the load 
was further increased gradually, a damage stabilization stage was entered when the 
additional damage growth was limited to the accumulation of damage in the previously 
experienced modes. At 8.99 kN (2,020 lb), two distinct laminate compressive fracture 
zones were formed on both sides of the loading nose in the top flange. In the bottom 
flange, 0° ply fiber tensile fractures grew into the five outermost 0° plies under the web, 
and grew into the four outermost 0° plies at nodes not located directly under the web. 
Fiber tensile fractures in the flange plies directly under the web were extended for 1 1 
nodes along the beam axis, centered at midspan. Approximately half the clear span length 
of the beam was affected by tensile fiber fractures in the center of the flange directly under 
the web. When the load was increased further, compressive laminate cracking began to 
grow rapidly in the upper flange, and fiber tensile fractures propagated along the beam 
axis in the bottom flange under the web, causing the total collapse of the beam. Figure 5.4 
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shows the beam finite element model immediately before collapse when some of the finite 
elements of the top and bottom flanges have been deleted due to the compressive 
fracturing of the top flange and tensile fracturing of the bottom flange. The simulated 
ultimate structural collapse mode was consistent with experimental observations reported 
in Takayanagi et al. (1994). However, the ultimate collapse mode is merely the final result 
of complex damage initiation, growth, and propagation characteristics that are determined 
via computational simulation. The fracture-resistant design of a laminated fiber composite 
beam requires a thorough evaluation of the damage progression characteristics, as well as 
the ultimate collapse load and the collapse mode. 

Figure 5.5 shows the relationship between load and deflection. The deflection coordinate 
plotted as the abscissa in Figure 5.5 represents the relative deflection of point A with 
respect to point B (see Figure 5.2). Experimental data, as well as the modified beam 
theory and composite beam theory results, shown in Figure 5.5 are from Takayanagi et al. 
(1994). Simulated displacements via the present method follow composite beam theory 
results very closely up to a load level of approximately 7 kN. It is noteworthy that the 
extensive matrix cracking and limited ply compressive failures that occurred prior to 
reaching the 7 kN load level did not affect the apparent structural response significantly 
under the monotonically increasing static loading. At 7.65 kN, when the outermost eight 
0° plies of the compression flange fractured, a substantial change occurred in the structural 
response that was observable in the load-deflection relationship. As the peak loading was 
approached, the simulated load-deflection relationship became more nonlinear compared 
to the experimentally observed load-deflection curve reported by Takayanagi et al. (1994). 
This difference may be explained as the effect of additional constraints imposed on the 
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upper flange by the loading yoke, which did not allow the free deformation of the 
fractured top flange when the loading was increased during testing. Nevertheless, the 
simulated peak load agrees well with that observed experimentally. 

In general, overall structural damage includes ply damage and also through-the-thickness 
fracture of the composite laminate. Varied and complex composite damage mechanisms 
are simulated via evaluation of the individual ply failure modes and associated degradation 
of laminate properties. The type of damage growth and the sequence of damage 
progression depend on the composite structure, loading, and material properties. A scalar 
damage variable, derived from the total volume of the composite material affected by the 
various damage mechanisms, is also evaluated as an indicator of the level of overall 
damage induced by loading. This scalar damage variable is useful for assessing the overall 
degradation of a given structure under a prescribed loading condition. The rate of increase 
in the overall damage during composite degradation may be used as a measure of 
structural propensity for fracture. The procedure by which the overall damage variable is 
computed is given by Minnetyan et al. (1990). 

Figure 5.6 shows the relationship between load and the damage energy release rate 
(DERR). The DERR is defined globally as the work done by applied forces, per unit 
damage produced during structural degradation. The DERR can be used to evaluate 
structural resistance against damage propagation at different stages of loading. A high 
DERR level indicates the dominance of global structural characteristics in damage 
progression. On the other hand, a low DERR level indicates the importance of local 
material properties on damage progression. The DERR usually starts at a high level at 
damage initiation. It is gradually lowered during damage growth and accumulation. Near 
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the peak load, DERR rates begin to increase again as the damage propagation stage is 
entered. The minimum DERR level signals the completion of the damage tolerant 
structural response stage. From Figure 5.6, the 8.10 kN load corresponding to the 
minimum value of the DERR indicates that the maximum load for damage tolerance was 
8.10 kN (1,820 lb). At 8.10 kN loading, there was some damage stabilization after the 
formation of through-the-thickness cracks in the top flange. Subsequent to the 8.10 kN 
load corresponding to minimum DERR, damage propagated more rapidly as the peak load 
was reached. 


Table 5.8 Relative Deflections under 4KN Loading 


Specimen 

Wobs (mm) 

Wmbt (mm) 

W C bt (mm) 

Wcod (mm) 

A 

0.384 

0.359 

0.348 

0.346 

B 

0.616 

0.621 

0.613 

0.631 

C 

1.136 

1.189 

1.184 

1.163 


Table 5.9 Peak Loads of CFRP Beams 


Specimen 

Failure 

Mode 

Peak Load (KN) 

Simulation 

Experiment 

A 

Compression flange 

9.54 

9.20 

B 

Upper web 

13.11 

12.78 

C 

Lateral bucking 

11.41 

11.82 


Figure 5.7 shows the relationship between load and damage. Damage started at 2.56 kN 
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(575 lb) loading. The amount of damage increased greatly after the load level of 8.10 kN 
was reached. When the load reached 9.52 kN (2,140 lb), the simulated specimen collapsed 
as the damage level reached 13.02%. The corresponding experimental collapse load was 
9.20 kN (2,068 lb). The experimentally observed peak load was about 3% lower than the 
simulated peak load. 

For specimen B, damage was initiated at the midspan of the top and bottom flanges 
simultaneously due to matrix tensile failures by the ply transverse tensile failure mode. The 
damage initiation load was 3.42 kN (776 lb). Damage propagated to 99 nodes as the 
loading increased to 3.97 kN (893 lb). When the load reached 12.32 kN (2,770 lb), most 
nodes were damaged and a through-the-thickness crack occurred at the upper middle part 
of the web due to ply compressive failures, as the reduced number of 0° flange plies 
caused a larger share of the load to be carried by the web rather than by the flanges. 
Subsequently, damage propagated to the middle of the upper flange and caused the I-beam 
to break completely. Figure 5.8 shows the simulated and experimentally observed load- 
deflection relationships for specimen B. From Figure 5.8, it may be noted that for 
specimen B the simulated relative deflections between points A and B are in very good 
agreement with the experimental results. The improved agreement of the simulated results 
with experimental data for specimen B may be explained by the fact that the quasi- 
isotropic outside flange laminate of specimen B is less critical than the web, and therefore 
is less sensitive to the loading and boundary conditions compared to the unidirectional 
laminate at the outer flange of specimen A. 

Due to the lack of 0° plies in the flanges, specimen C was susceptible to large vertical 
deflections and lateral buckling. A finite element structural stability analysis for specimen 
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C indicated that lateral buckling would occur under a 11.41 kN (2,565 lb) load. The 
experimentally observed buckling load for specimen C was 11.82 kN (2,658 lb). 
Computational simulation indicated that stress-induced composite damage did not occur 
prior to reaching the buckling of the beam. Specimen C was considerably more flexible 
compared to specimens A and B. The relative deflections between points A and B (see 
Figure 5.2), for specimens A, B, and C, under a load of 4 kN, obtained via computational 
simulation (W MB t), modified beam theory (W CB t,)> composite beam theory (W CB t,X and 
experimentally observed data (W obs ), are shown in Table 5.8. 

Table 5.9 shows the peak loads obtained by computational simulation and experimental 
data from Takayanagi et al. (1994) for the three specimens. For all three specimens, the 
maximum difference between the simulated peak load and corresponding experimental 
failure observations was less than 3.5%. From the simulated and observed fracture modes, 
it can be concluded that for specimen A the peak load may be increased by adding 90° 
plies to the top flange to prevent the transverse tensile fracture of the 0° plies. For 
specimen B, additional 0° plies in the flange would reduce the stress concentration in the 
web and thus increase the peak load. For specimen C it would be necessary to add 0° plies 
into the flanges to prevent lateral buckling. 
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Figure 5.4 Finite Element Model of Specimen A at Collapse 


5.6 Conclusions 

In light of the damage propagation evaluation of the CFRP I-beams and from the general 
perspective of the integrated computational evaluation method, the following conclusions 
are drawn: 

1. The assessment and design of CFRP structural shapes may be improved by an 
evaluation of the composite damage initiation and damage propagation mechanisms. 

2. In general, composite structures demand a detailed understanding of the stress states 
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and material degradation characteristics for effective design. 

3. Computational simulation, with the use of established composite mechanics and finite 
element methods, can be used to simulate the various stages of degradation in a composite 
structure, such as damage initiation, damage growth, damage accumulation, damage 
propagation, and structural fracture. 

4. Computational simulation is based on constituent level material properties. However, 
even if the fiber and matrix properties are not explicitly defined, they can be determined 
from the given composite ply properties using a least-squares minimization procedure. 

5. Significant design parameters such as composite constituent stiffness, strength, and 
effective ply configurations may be identified with the help of a composite mechanics 
module and computational simulation. 

6. The demonstrated procedure is flexible and applicable to all types of constituent 
materials, structural geometry, and loading. Hybrid composites and homogeneous 
materials, as well as binary composites can be simulated. 

7. The present computational simulation procedure provides a new general methodology 
to investigate damage propagation and progressive fracture for any structure. 
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Figure 5.5 Relative Deflection of Specimen A under Loading 
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Figure 5.6 Damage Energy Release Rates of Specimen A under Loading 
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Figure 5.7 Damage Progression of Specimen A under Loading 



Figure 5.8 Relative Deflection of Specimen B Under Loading 
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Chapter 6 


Computational Simulation of 
Progressive Fracture in Stitched 

J-Stiffened Composite Shear Panels 
in the Postbuckling Range 

6.1 Introduction 

Advanced composite structures, such as three-dimensional woven and stitched stiffened 
composite laminates are used in many aerospace applications such as aircraft fuselage and 
wing structures. Stitching or other through the thickness reinforcement greatly reduces the 
driving force for propagation of the delamination crack. The stitching creates closing 
traction acting across the laminate thickness that limits damage extension. However, it is 
important that the proportion of the through the thickness reinforcement is minimized to 
avoid significant degradation of the in-plane stiffness and strength properties. 

The design evaluation of composite structures requires an assessment of their safety and 
durability under service loads and possible overload conditions. However, the architecture 
of a three-dimensionally reinforced composite is complex and, therefore, the parameters 
controlling its mechanical properties are numerous. Mechanical testing to characterize all 
the effects of composite architecture could be an economically unrealistic proposition. It is 
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therefore crucial to develop experimentally verified analytical models to predict the load 
carrying capacity, damage tolerance, safety, and reliability of composite structures. An 
integrated computer code was developed in this paper for the purpose of evaluating 
stiffened composite panel progressive fracture in the post-buckled range. 

An important feature of computational simulation is the assessment of damage stability or 
damage tolerance of a structure under loading. At any stage of damage progression, if 
there is a high level of structural resistance to damage progression under the service 
loading, the structure is stable with regard to fracture. The corresponding state of 
structural damage is referred to as stable damage. On the other hand, if damage 
progression does not encounter significant structural resistance, it corresponds to an 
unstable damage state. Unstable damage progression is characterized by very large 
increases in the amount of damage due to small increases in loading, whereas during stable 
damage progression the amount of increase in damage is consistent with the increase in 
loading. 

Internal damage in composites is often initiated as cracking due to normal stress transverse 
to fiber orientation. Further degradation is in the form of additional cracking, 
delamination, and fiber fractures that may lead to structural fracture. The novelty of this 
paper is that various stages of degradation, from damage initiation to structural collapse of 
stitched laminates are simulated for the first time. The present computational method has 
previously been used for the investigation of the behavior of conventional laminated 
composite structures, including the effects of composite degradation on structural 
response (Minnetyan et al 1991), composite structures global fracture toughness 
(Minnetyan et al 1990), effect of the hygrothermal environment on durability (Minnetyan 
et al 1992-1), damage progression in composite shells subjected to internal pressure 
(Minnetyan et al 1992-2), the behavior of discontinuously stiffened composite panels 
under loading (Gotsis et al 1996), damage propagation in thick composite shells under 
external pressure (Minnetyan and Chamis 1997), damage and fracture progression in 
CFRP I-beams with three different lay-ups (Huang and Minnetyan 1998). 


NASA/CR— 2005-213963 


82 



6.2 Methodology 

The computational simulation tool presented in this paper is an integrated computer code 
consisting of three modules: composite mechanics, finite element analysis, and damage 
progression modeling. The overall evaluation of composite structural durability is carried 
out in the damage progression module that keeps track of composite degradation for the 
entire structure. The damage progression module relies on the composite mechanics 
module for micromechanics and macromechanics, and calls a finite element analysis 
module to obtain generalized stresses and displacements for the composite structure. 

6.2.1 Evaluation of Elastic Constants 

Elastic constants are essential for stress analysis of composites. The stiffness averaging 
method that was developed by Kregers and Melbardis (1978) is widely used to predict the 
deformation characteristics of a composite with three-dimensional reinforcement from the 
known mechanical properties of its components. However, utility of the stiffness 
averaging method is limited to the prediction of overall stiffness properties of a structure 
for elastic analysis. Stiffness averaging does not retain the information on the spatial 
configurations of fiber reinforcements. Therefore, the effects of damage and degradation 
cannot be directly taken into account for a specific ply and the results of structural analysis 
cannot be decomposed to the micromechanics level for detailed analysis of damage 
progression. In this paper, the stiffness averaging method and classical laminate theory are 
combined together to predict the properties for 3-D reinforced composites. For stitched 
laminates, each ply consists of not only in-plane fibers, but also through-the-thickness 
stitching fibers. Accordingly, each ply is divided into two subvolumes. The first 
subvolume consists of only in-plane fibers and the second consists of only the stitching 
fibers. The properties of each subvolume in its local coordinate system are computed from 
the properties of fiber and matrix, based on the classic intraply composite micromechanics 
theory (Chamis and Sinclair 1979), and then transformed to structural directions. The ply 
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properties are obtained by adding contribution of each subvolume based on the stiffness 
averaging method: 



Where Ec is the composite ply stiffness. Vi is the calculated volume of the i-th subvolume, 
N is the number of subvolumes. E, is the stiffness for the i-th subvolume, i?, is the 
coordinate transformation matrix for the i-th subvolume. After computation of the layer 
properties from the properties of subvolumes, properties of the entire stitched laminate are 
obtained by using laminate theory as one would obtain laminate properties for unstitched 
composites. 

6.2.2 Finite Element Analysis 

After the micromechanics analysis module generates the elastic properties for the stitched 
composite, a finite element analysis module is called to analyze the entire structural 
response. In general, the type of finite element model used depends on the complexity of 
the structure and the availability of computer resources. There are two possible choices 
for the analysis of composite structures. One is using anisotropic three-dimensional solid 
elements such as hexahedral or brick elements that accept the computed three dimensional 
composite properties directly. However, the modeling of a practical composite structure 
with three-dimensional elements requires huge computer resources and it is usually 
impractical for analyzing real structures. The second option is to use anisotropic shell 
elements that use the composite plate/shell element properties. Anisotropic plate or shell 
elements represent through-the-thickness properties of the stitched composite and are 
more efficient than three-dimensional elements. Therefore, implementation in this paper is 
focussed on the use of plate/shell elements in the finite element analysis module. The finite 
element module accepts the force-deformation relations computed by the composite 
macromechanics module, and carries out a stress analysis to generate the generalized 
stresses Nx, Ny, Nxy, Mx, My, Mxy, Qxz, Qyz for each node. The generalized stresses are 


NASA/CR— 2005-213963 


84 



supplied back to the composite mechanics module for the computation of local ply/stitch 
stresses and failure analysis. 


6.2.3 Failure Criteria 

Progressive damage and fracture simulations are carried out by imposing failure criteria 
locally within each subvolume with reference to the local coordinate orientations in the 
material directions. At each individual load step, stresses in stitching and in-plane 
subvolumes that can be obtained through the composite microstress analysis are assessed 
according to distinct failure criteria. The first twelve modes of failure are associated with 
the positive and negative limits of the six local stress components in the material directions 
as follows: 

^ (l 1C < °fll < S euT 
S (22 C ^ ^ (22 < ^ (22T 
S(33C < °/33 < S (3 3T 

S(12(-) < < ^12(+) 

S(23(~) < ®(23 < ^(2X + ) 

^ 03 (-) < ®(\3 < ^ 13 (+) 

The thirteenth failure mode is a combined stress or modified distortion energy (MDE) 
failure criterion that is obtained by modifying the usual distortion energy failure criterion. 
The modification takes into account the significant differences in the stress limits of the 
longitudinal and transverse directions of an orthotropic composite ply. Each component 
of ply stress is normalized with respect to its limiting strength. No relationship is assumed 
between normal and shear strengths. The MDE failure criterion has been demonstrated to 
be a good predictor of combined stress failure in composites. It may be considered as a 
variation of the Tsai-Hill theory (Tsai 1968). The MDE failure criterion (Chamis 1969) 
can be expressed as: 
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Where a and /f indicate tensile or compressive stresses, S eUa is the local longitudinal 
strength in tension or compression, S tl2a is the transverse strength <in tension or 
compression, and the directional interaction factor K ma p is defined as: 

!* + ~~ ^13 ) L , 22 0 ; ^23 ’ J 

n ^~[E n E n ( 2+v n +v„Y2 + v I , 

The directional interaction factor reduces to unity for homogeneous isotropic materials. 


6.2.4 Simulation of Damage Progression 

After each finite element stress analysis, failure criteria are used to evaluate possible failure 
within each subvolume of each ply at each node of the composite structure. Once the 
damage modes at each node are assessed, a damage index is created to record the damage 
information for each damaged node. The damage index contains the node number, the ply 
number, and the list of damage criteria that have become activated. When a new failure 
occurs at a subvolume, the damage index is updated correspondingly. The properties of 
each domain are degraded according to their damage index. 

If there is no damage after a load increment, the structure is considered to be in 
equilibrium and an additional load increment is applied. If any damage occurs, the 
composite properties affected by the damage are degraded, the computational model is 
reconstituted with updated finite element mesh and material properties, and the structure is 
reanalyzed under the same load increment. After reanalyzing, if there is any additional 
damage, the properties are degraded again and the structure is reanalyzed again. This 
cycle continues until no further damage occurs. 

The damage progression module keeps a detailed account of composite degradation for 
the entire structure. It also acts as the master executive module that directs the composite 
mechanics module to perform micromechanics and macromechanics analysis/synthesis 
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functions, and calls the finite element module with thick shell analysis capability to model 
stitched laminated composites for global structural response. 


6.2.5 Damage Energy Release Rate 

The measure of global fracture toughness is defined in terms of the Total Energy Release 
Rate (TERR) that is equal to the amount of damage energy expanded for the creation of 
unit damage volume in the composite structure. The magnitude of TERR varies during 
progressive degradation of the composite structure under loading, reflecting the changes 
in the fracture toughness of the composite. Computation of TERR during progressive 
fracture is needed in order to evaluate the composite fracture toughness and the degree of 
imminence of failure. 

TERR is computed as the ratio of damage energy to the corresponding damage volume 
that is generated: 

TERR=Damage Energy/Damage Volume 

Assuming locally linear elastic properties prior to damage, the damage energy can be 
computed as follows: 

Damage Energy — m 0.5S]/E j )v i 
‘ j 

Where Sj is the local composite strength associated with the damage mode, £} is the elastic 
modulus corresponding to Sj, and V, is the volume of damage. Typically, at the stage of 
damage initiation, there is a high rate of energy release that dissipates a significant portion 
of the strain energy stored in the composite structure. After the first burst of energy 
release, TERR usually drops down to a lower level, indicating the significant reduction in 
fracture stability of the damaged composite under the applied loading. TERR usually 
reaches its peak value as global structural failure occurs. 
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6.3 J-Stiffened Composite Panels 

The panels simulated in this paper contain two J-shaped stiffeners as shown in Figured. 1 . 
The experimental results for these panels were adopted from the literature (Yeh and Chen 
1996). The panels were made of AS-4/3501-6 carbon/epoxy unidirectional fabric. The 
fabric consisted of a wide sheet of unidirectional tows of fibers based together with 
polyester thread to keep it from unraveling. Twelve plies of the skin were laid up to form a 
1.83 mm (0.072 in.) thick quasi-isotropic [0,90,45,0,-45 ,90] s laminate, each ply had 
approximately the same thickness. The “J” shaped stiffeners were constructed with the 
same lay-up except for the flanges that were only half the thickness of the basic laminate. 
The skin panel and the stiffeners were stitched together with Kevlar and fiberglass threads. 

Yeh & Chen (1996) tested two panel specimens with aluminum doublers bonded onto 
edges. Panel #1 is shown in Figure 6.1(a) and Panel #2 is shown in Figure 6.1(b). In order 
to simulate these specimens, a finite element model shown in Figure 6.2 was used. Mindlin 
type thick composite shell elements were utilized to represent the skin panel and stiffeners. 
The aluminum doublers were modeled with rigid elements. 


6.4 Fiber and Matrix Properties Calibration 

The first step in this simulation was to identify details of the fiber and matrix properties. 
Since Yeh & Chen (1996) only gave ply properties, and did not provide explicit 
information on fiber and matrix properties, the composite mechanics module (Murthy and 
Chamis 1986) was used to calibrate the fiber and matrix properties so that the computed 
ply properties matched those of the AS-4/3501-6 composite tested by Yeh & Chen. The 
properties for AS-4 and 3501 provided by the resident databank of the composite 
mechanics code (Murthy and Chamis 1986) were used as the initial values of AS-4 fiber 
and 3501-6 matrix respectively. Ply properties computed by the composite mechanics 
module were compared with the experimental monolayer values provided by Yeh & Chen. 
The difference between the computed properties and experimental observation was 
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Figure 6.2. Finite Element Model and Loading Condition 



Figure 6.3. Panel Buckling Mode 
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Where p i ex P enmen ' i s the value of composite property (such as Eu, E 22 , Gn , etc.) from 
experiment. p. com P uted j s the value of the same property computed by the composite 
mechanics module. N is the number of composite ply properties observed. 

It is obvious that the smaller the value of Q the more accurate the fiber and matrix 
properties. However, because of the natural complexity of the composite material, as well 
as the uncertainty and variability of experimental observations, it was not possible to 
obtain computed ply properties exactly the same as those experimentally observed. Based 
on the fiber and matrix properties given in Tables 6.1 and 6.2, the composite properties 
computed by the composite mechanics module are shown in Tables 6.3 and 6.4. The 
value of Q for the monolayer properties is 0.0272. This means that the overall difference 
between experimental and simulated ply properties is less than 3 percent. 


Table 6. 1 : Fiber Properties 



AS4 

Kevlar 

Number of fibers per end 

10000 

580 

Fiber diameter (in) 

0.0003 

0.00046 

Fiber Density ( lb/in) 

0.063 

0.053 

Longitudinal normal modulus (l(f psi) 

31.2 

22.0 

Transverse normal modulus (10 6 psi) 

2.82 

0.60 

Poisson's ratio |Xi 2 

0.27 

0.35 

Poisson's ratio |i 2 3 

0.30 

0.35 

Shear modulus Gi 2 (10 6 psi) 

10.0 

0.42 

Shear modulus G 23 (10 6 psi) 

10.0 

0.22 

Longitudinal thermal expansion coefficient (10' 6 /°F) 

0.550 

2.20 

Transverse thermal expansion coefficient (10' 6 /°F) 

0.560 

30.0 

Longitudinal heat conductivity ( BTU- in/hr/in /°F) 

4.03 

1.70 

Transverse heat conductivity ( BTU -in/hr/in /°F ) 

4.03 

1.70 

Heat capacity ( BTU/lb/°F) 

0.17^ 

0.25 

Tensile strength (ksi) 

320 

400 

Compressive strength (ksi) 

295 

75 


NASA/CR— 2005-213963 


91 






Table 6.2. 3501-6 Matrix properties 


Matrix density (Ib/in 1 

0.0457 

Normal modulus (ksi) 

600 

Poisson's ratio 

0.400 

Coefficient of thermal expansion (10' 4 /°F) 

0.4 

Heat conductivity ( 10' 3 BTU -in/hr/in /°F ) 

8.681 

Heat capacity (BTU/lb/°F) 

0.250 

Tensile strength ( ksi) 

4.2 

Compressive strength ( ksi) 

3.68 

Shear strength ( ksi) 

14.0 

Allowable tensile strain 

0.02 

Allowable compressive strain 

0.05 

Allowable shear strain 

0.04 

Allowable torsional strain 

0.04 

Void conductivity ( BTU -in/hr/in 2 /°F) 

0.225 

Glass transition temperature (°F) 

420 


Table 6.3. Monolayer properties and allowables for AS4/3501-6 unidirectional fabric 


Experiment 

Simulation 

longitudinal modulus El (psi ) 

17.5 x 10 6 

17.43 x 10 6 

transverse modulus Er(psi) 

1.45 x 106 

1.42 x 106 

shear modulus Glt (psi) 

0.51 x 106 

0.51 x 106 

poisson’s ratio v L t 

0.34 

0.33 

longitudinal tension strength Su (psi) 

176,000 

176,000 

longitudinal compression strength Su (psi) 

162,750 

162,200 

transverse tension strength Sn (psi) 

3,493 

3,566 

transverse compression strength Sn (psi) 

31,388 

31,250 

shear strength S s h (psi) 

11,547 

11,890 


6.5 Results and Discussion 

Figure 6.2 shows the finite element model used in the simulation of these panels. The 
model is loaded by applied tension loading in the diagonal direction at one comer of the 
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panel, and constrained in X, Y, and Z directions at the opposite corner. First, the 
simulation under a monotonically increasing static loading was performed, without 
considering lateral buckling of the skin. Since experimental data indicated that local 
buckling of the skin occurred prior to reaching the peak load, a linear buckling analysis 
was conducted next. Buckling analysis showed that the initial buckling load was 8275 lb. 
That is very close to the experimental results. According to Yeh and Chen’s experiments 
(1996), initial buckling load for panel #1 was approximately 75001b and for panel #2 was 
approximately 8500 lb. In order to study the effect of lay-up configuration on the peak 
load, four different lay-ups described in Table 6.4 were simulated in the present study. 
Since Yeh and Chen gave experimental results for the [0/90/45/0/-45/90] layup, the 
simulated details of damage progression for this lay-up are presented in the next section. 
Only a summary of the simulation results will be presented for panels with other lay-ups. 


6.5.1 Analysis Without Considering Lateral Buckling of the Skin 

In order to study how buckling can affect damage progression, first the stiffened panel 
with [0/90/45/0/-45/90] s lay-up was analyzed without considering the lateral buckling 
effect. The load applied to the panel was started at 445 N (100 lbs) and then increased 
gradually. When the load increased to 133 kN (29900 lbs), damage occurred at the edges 
of the stiffener due to the transverse compression failure of some 0-degree plies (plies 
9,12,13,16). These damaged nodes are indicated by an “a” in Figure 6.4. When the load 
increased to 139 kN (31300 lbs), damage occurred at another end of the edge of the 
stiffeners, indicated as “b” nodes in Figure 6.4. The “b” nodes were also damaged due to 
compression failures in the 0-degree plies. When the load was increased to 152 kN (34200 
lbs), the damaged nodes extended to the middle from the ends along the edges of the 
stiffeners and 90-degree plies started to damage due to shear failure, as indicated by the 
“c’ nodes in Figure 6.4. When the load reached 165 kN (37000 lbs), there were more 
nodes damaged, indicated as “d” nodes in Figure 6.4. 165 kN was the maximum 
equilibrium load before the panel collapsed. When load reached 169 kN (37900 lbs), 
some nodes shown in Figure 6.4 as “x” nodes started to fracture and nodal fractures 
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immediately extended to adjacent nodes along the edges of the J-stiffeners. Even though 
the load was not increased any more, the fracture line continued to extend causing the 
collapse of the entire structure. 


Table 6.4. Laminate properties for panels with different lay-up 



[0/90/45/0/-45/90] s 
Experiment simulation 

[45/-45/90 

/45/0/-45] s 

[0] 12 

[0/90/0/90 

/0/90] s 

E x (10 6 psi) 

7.765 

7.78 

5.40 

17.43 

9.497 

E y ( 10 6 psi) 

7.272 

7.78 

5.40 

1.42 

9.497 

Gxy( 10 6 psi) 

1.95 

1.85 

3.186 

0.51 

0.51 


Table 6.5. Initial damage and peak loads for panel with different lay-up 


Lay-up 

Initial 

damage 

load 

(lb) 

Maximum 
equilibrium 
load (lb) 

Load causing 
structure 
collapse(lb) 

Average 
Experimental 
failure load(lb) 

[0/90/45/0/-45/90] s 

12800 

27600 

27800 

27500 

[45/-45/90/45/0/-45] s 

9685 

26000 

26300 

- 

[0/90] 6s 

3780 

8740 

8810 

- 

[0] 12 

850 

6570 

6600 

- 


6.5.2 Postbuckling Analysis 

Since the initial local buckling load is much less than the peak load, post-buckling analysis 
is necessary for a realistic simulation of the stiffened panel. In order to do postbuckling 
analysis, buckling analysis was conducted before progressive fracture to obtain the 
buckling load. Subsequently, damage progression was simulated as follows: If the applied 
load was less than the buckling load, then a regular FEM analysis was conducted for the 
evaluation of composite stresses. After the applied load exceeded the buckling load, the 
buckling mode as shown in Figure 6.3 was imposed onto the structural geometry for FEM 
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analysis. The load was started at 445 N (100 lbs). When the load increased to 36.8 kN 
(8275 lbs), local buckling occurred, therefore the buckling mode as shown in Figure 6.3 
was used for the FEM analysis increments afterwards. When load was increased to 56.9 
kN (12800 lbs), the “a” node shown in Figure 6.5 started to sustain damage due to the 
transverse compression failure of the outmost 0-degree ply. When the load increased to 
66.7 kN (15000 lbs), damage extended to the edges of the stiffeners and other nodes (see 
“b” nodes in Figure 6.5). As the load was increased, damage extended to more nodes. 
When the load increased to 123 kN (27600 lbs), 90-degree plies were damaged in the 
shear failure mode. However, no node had yet experienced a through-the-thickness 
laminate fracture. The positions of nodes damaged at this stage are indicated as “c” nodes 
in Figure 6.5. When the load was increased to 124 kN (27800 lbs), through-the-thickness 
laminate fracture occurred at some nodes indicated as “x” nodes in Figure 6.5. Through- 
the-thickness fracture extended to adjacent nodes very quickly, resulting in structural 
fracture. Therefore the peak load was 27600 lbs. This was very close to the experimental 
peak load. According to Yeh & Chen (1996), the experimental peak load for Panel #1 was 
119 kN (26800 lbs), and for Panel #2 was 125 kN (28200 lbs). The average experimental 
peak load was 122 kN (27500 lbs). Experimental and simulated relationships between the 
applied load and the corresponding diagonal deflection are shown in Figure 6.6. 

Figure 6.7 shows the relationship between applied load and produced damage volume. 
Damage initiation occurred at 56.9 kN (12800 lbs). After the load reached 109 kN 
(24500 lbs), the amount of damage increased greatly. When the load reached the 123 kN 
(27600 lbs), 5.84 percent of the panel was damaged. After that, the panel entered an 
unstable damage propagation state. When the simulated load was increased to 124 kN 
(27800 lbs), the damage volume increased rapidly until peak structural failure occurred. 
Figure 6.8 shows the relationship between damage energy and load. Similar to the damage 
volume depicted in Figure 6.7, damage energy also increased rapidly after the load reached 
109 kN (24500 lbs). 
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6.5.3 Analysis of Panels with other lay-ups 

In order to study how the laminate configuration affects the peak load, simulations of 
panels with other lay-ups were also included in this study. All the panels had the same size, 
same thickness and same stiffeners. The three additional simulated panels had the laminate 
configurations of [45/-45/90/45/0/-45] s , [0/90] 6s , [0] i 2 - The laminate properties for these 
panels are shown in Table 6.4. Tables 6.4 and 6.5 show the initial damage load, peak load, 
and initial buckling load. From these results it can be seen that the lay-up configuration has 
a significant effect on the peak load. Similar to the panel with a lay-up [0/90/45/0/-45/90] s , 
the panel with a lay-up of [45/-45/90/45/0/-45] s was locally buckled before its collapse. 
After local buckling, the load was able to increase quite a bit. The peak load was more 
than three times the initial buckling load. For the panels with a cross-ply lay-up of [0/90] 6s , 
the peak load was just a little larger than the initial buckling load. This should not be 
interpreted to mean that buckling has a much more significant effect on the peak load for 
the cross ply panel than for the panels with 45-degree plies. As a matter of fact, even when 
the initial buckling effect was not taken into account, the peak load for the cross ply panel 
was only 42.1 kN (9470 lbs). Therefore the 45-degree plies play an important role in 
increasing peak loads. For the panel with unidirectional lay-up, the initial buckling load 
was larger than the peak load, therefore the unidirectional panel collapsed before buckling. 
Figure 6.9 shows the load-deflection relationships for the panels with different ply layup 
configurations. 


6.6 Conclusions 

1. A Computational simulation tool developed for the analysis of conventional laminated 
composite structures (Minnetyan et al, 1992, 1995) was extended to the analysis of 
composite-stiffened-stitched-shear panels. 

2. Postbuckling effects have a significant influence on the peak collapse load of 
quasiisotropic stiffened shear panels. 
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3. Computational simulation without considering the postbuckling mode overpredicts the 
collapse load by approximately 35 percent. The stiffness is overpredicted by 7.6 percent. 

4. It was found that ply configuration significantly affects the strength of stitched 
composite panels under shear loading. 

5. Damage volume and/or damage energy can be used to evaluate the overall degradation 
stage of a structure. 

6. Computational simulation provides an effective means for parametric design 
investigations of stiffened and stitched composite structures. 



Note: Initial damage load for nodes a: 128001b, b: 150001b, c: 276001b, 
x indicates the nodes first cracked under loading 379001b. 

Figure 6.4. Damage Node Indication for Panel with Ply Configuration [0/90/45/0/-45/90]s 
with Postbuckling Analysis 
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Note: Initial damage load for nodes a: 299001b, b: 313001b, c: 342001b, d: 370001b, 
x indicated the nodes first cracked under loading 379001b 


Figure 6.5. Damage Node Indication for Panel with Ply Configuration [0/90/45/0/-45/90] s 
without Postbuckling Analysis 
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Chapter 7 


Computational Simulation of 
Progressive Fracture in 
Ceramic Matrix Composites 

7.1 Introduction 

Ceramic matrix composites (CMC) are currently a subject of great deal of research 
interest due to the requirement of new high temperature materials for use in high-speed 
engine structural components and other applications. Continuous-fiber reinforced ceramic 
matrix composites offer a great deal of promise in this area and hence considerable effort 
is being devoted for their development. They offer several advantages such as high specific 
stiffness and strength, high toughness as compared to monolithic ceramics, environmental 
stability and wear resistance for both room and elevated temperature applications. Ceramic 
matrix composites are reinforced primarily for enhancing the toughness as the matrices in 
these composites are brittle and fail at relatively low strain levels. It is also desirable to 
have a weak fiber/matrix interface in these composites to enable such toughening 
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mechanisms as fiber debonding, fiber bridging, fiber pullout, crack deflection etc. to 
eliminate catastrophic failures. However, for the lull potential of these materials to be 
attained, research will have to provide insights into accurate micro-mechanical 
characterization of ceramic matrix composites. This chapter will develop computational 
tools to simulate the progressive fracture in ceramic matrix composite structures. 

7.2 Micro- and Macro-mechanical Characterization of 
Ceramic Matrix Composites 

Prediction of properties of ceramic matrix composites is based on a computational method 
developed to specifically simulate aspects unique to ceramic matrix composites (Mital and 
Murthy 1996). The method incorporates a slice-by-slice substructuring technique that has 
many advantages including a more accurate micro-mechanical representation of interfacial 
conditions, both around the fiber and through-the-thickness and provides a greater detail 
of stress distribution within a ply. 

A multilevel substructuring technique, which includes a unique fiber substructuring 
concept, is used for the analysis of continuous fiber reinforced ceramic matrix composites. 
This technique has four levels of substructuring-from laminate to ply, to subply, and then 
to fiber. The fiber is substructured into several slices and the micromechanics equations 
are applied at the slice level. Although the basic philosophy can be applied to the analysis 
of any continuous fiber reinforced composite, the emphasis is on the development of a 
method to specifically analyze and simulate aspects unique to ceramic matrix composites. 
The aspects of interest include varying degrees of interfacial bond around the fiber 
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circumference and accounting for the fiber breaks and local matrix cracking which may 
lead to rapid degradation of interphase at higher temperatures due to oxidation. In 
addition, the multilevel substructuring technique can account for different fiber shapes and 
integrate the effect of all of these aspects on composite properties/ response and, in turn, 
provide a greater detail in stress distribution. 

Composite micromechanics equations are used to determine the equivalent elastic 
properties of a composite material in terms of the elastic properties of the fiber and matrix 
constituent materials. The properties of interest are composite moduli, Poisson's ratios, 
thermal expansion coefficients, thermal conductivity, heat capacity, etc., and various 
composite strengths. The micromechanics equations are derived for a representative 
volume element (RVE), sometimes referred to as a "Unit cell" which is the smallest region 
over which the stresses and strains are assumed to be macroscopically uniform. The unit 
cell consists of fiber, matrix, and possibly an interphase treated as a separate constituent. 
The geometry of the unit cell depends upon the chosen array pattern for the fibers, e.g., 
square, hexagonal or any other kind of repeating geometry. Equivalent properties for the 
unit cell or RVE are then derived based on the constituent properties using a mechanics of 
materials approach. This approach is based on certain assumptions, in addition to a chosen 
array pattern. The assumptions are: (a) all the constituents are subjected to the same 
strain in the fiber direction in the case of a unidirectional composite, and (b) the same 
transverse stress is applied to all constituents in the direction transverse to the fiber. These 
are standard assumptions used for the derivation of micromechanics equations, but they 
are not mathematically rigorous. A mathematically rigorous solution that ensures the 
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displacement continuity across the fiber and matrix boundary can be obtained through the 
use of theory of elasticity. 

Prior to fiber substructuring the unit cell may consist of three distinct regions A, B, and C 
(figure 7.1 A). The region A consists of pure matrix, the region B consists of matrix and 
interphase, and region C consists of fiber, matrix, and interphase. The interphase is treated 
as a distinct region with distinct mechanical and thermal properties. Thus, it can represent 
a reaction zone formed due to the chemical reaction between fiber and matrix or a 
separate layer provided intentionally to prevent such a reaction. The different regions 
facilitate representation of nonuniformity in the local stress distribution. 

This approach has been taken one step further for the analysis of ceramic matrix 
composites. Although the same square fiber array pattern is assumed, the unit cell is 
further subdivided into several slices. The equations of micromechanics are derived for 
slices, i.e., slice equivalent properties are computed based on the properties of the fiber, 
matrix, and interphase. The fiber substructuring and slice geometry are shown in figure 
7. IB. The derivations of the micromechanics are provided in detail by Murthy and Chamis 
(1992). For example, if k f , k m and h are fiber, matrix, and interphase volume ratios 
respectively, then for a typical slice 


U f l 

k =—■ k —1- 


V 

s s 

Where s is the width of unit cell, d f is the width of fiber, d m is the width of matrix, di is the 


k s =2— 
s 


width of interface. 

The equivalent longitudinal slice modulus is given by, 

E/ll = kf Efll ^ mil 
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A 

B 


C 


B 
A 

A). Unit cell square array concepts of ceramics composite 



B). Horizontal slicing 


C). Vertical slicing 


Figure 7. 1 Substructuring for ceramics composite micromechanics 
The transverse modulus in the 2-2 direction is given by. 
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and the longitudinal thermal expansion coefficient is given by, 

k f (Xfu E fu &mii E m ii + k t a jU E in 

a ‘" = E 

11 

Similarly, other mechanical and thermal properties are expressed in terms of the 
constituent properties. Once the equivalent slice properties are obtained, the equivalent 
properties of the unit cell or RVE are obtained by using laminate theory in an analogous 
manner as one would obtain laminate properties from ply properties. It should also be 
mentioned that 2-2 or horizontal slicing is used to compute 1-1 and 2-2 slice properties, 
while 3-3 or vertical slicing is used to compute slice properties in the 3-3 direction. 

To obtain the equivalent properties of a unit cell or RVE, the slice properties are 
integrated by using classical laminate theory. If there is only one fiber through the ply 
thickness, then the ply equivalent properties are identical to the properties of the unit cell 
or the RVE. If there are a number of fibers through the thickness of a single ply, then the 
unit cell properties have to be integrated by using the laminate theories again to obtain the 
ply properties. Laminate or composite properties are obtained from the individual ply 
properties by using macromechanics theories. For braided or woven composites, 
properties of the out-of-plane fiber reinforcements can be rotated into the laminate 
structural axes based on coordinate transformations and integrated into the entire structure 
via the zone-by-zone stiffness averaging method presented in chapter 3. 

After the laminate properties are computed, the next step is to obtain the laminate 
response due to externally applied loads. The flowchart to outline these procedures is 
shown in figure 7.2. The left side of the chart shows synthesis of the properties from slice 
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to laminate level. The right hand side of the flowchart shows the decomposition of the 
response from laminate to ply, to slice, and then to constituent microstress levels. 



Figure 7.2 Integrated analysis procedure for ceramic matrix composites 


7.3 Failure Criteria for Ceramic Matrix Composites 


Failure criteria are essential for predicting progressive fracture in composites. Ceramic 
matrix composites consist of fiber, matrix, and interface. Fracture initiation in composites 
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is usually caused by fiber breaking, matrix cracking, and/or interface failure. Therefore, in 
order to predict composite failure, we need to predict constituent failures. The failure 
criteria associated with stress limits can be expressed for fiber stresses as follows: 

S fuc < a f ii < ^ fin 

S f22C & f 22 ^ $ f22T 

S flic < G fV> < Sf33T 
S fl2S(,~) < 'C fl2 < $ fl2S(+) 

S f23S(-) < T /23 < *^/23S(+) 

*^/13S(— ) < T /13 < ^/135(+) 

Where 

Sfin — Longitudinal tensile strength of fiber, 

Sfuc — Longitudinal compressive strength of fiber, 

S/22T — Transverse tensile strength (Y direction) of fiber, 

Sf22c — Transverse compressive strength (Y direction) of fiber, 

Sf33T — Transverse tensile strength (Z direction) of fiber, 

Sf33c — Transverse compressive strength (Z direction) of fiber, 

Sfi2s — XY plane shear strength of fiber, 

Sfus — XZ plane shear strength of fiber, 

Sf23s — YZ plane shear strength of fiber 

Failure criteria for matrix and interface can be written similarly. 


NASA/CR— 2005-213963 


108 



7.4 Progressive Fracture Simulation 


An integrated computer code is developed to simulate progressive fracture in fiber 
reinforced ceramic matrix composite materials. As such it facilitates the evaluation of 
ceramic matrix composite structures with regard to damage tolerance and durability by 
simulating the various stages of degradation, from damage initiation to complete failure. 
This code contains three modules: CMC (Ceramic Matrix Composite micro-mechanics 
and macromechanics) module, FEM (Finite Element Method) module, and PFA 


(Progressive Fracture Analysis) module. Figure 7-3 shows schematic flowchart of this 
code. 


Damage quantification 




Figure 7.3 Ceramic Matrix Composite Progressive Simulation Cycle 
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The CMC module is based on software developed to specifically simulate aspects unique 
to ceramic matrix composites (Mital and Murthy 1996). It incorporates a fiber 
substructuring technique that provides: 1) accurate micromechanical representation of 
interfacial conditions, both around the fiber and through-the-thickness, and 2) greater 
detail of stress distribution within a ply. This module generates the properties (A, B, D 
matrix) of ceramic matrix composites from properties of the constituents (fiber, matrix, 
and fiber-matrix interface). The generated structural properties are used in the FEM 
module. CMC module can also be used to determine constituent level stresses, such as 
fiber, matrix and interface stresses. 

The FEM module performs finite element analysis to determine the generalized stresses 
over the entire ceramic matrix composite structure. These generalized stresses are referred 
back to the CMC module for iterative determinations of composite microstresses. 

The PFA module keeps a detailed account of composite degradation over the entire 
structure and also acts as an executive control module to direct the CMC module to 
perform micromechanics, macromechanics analysis, and synthesis functions. The damage 
progression module also calls FEM module with anisotropic thick shell analysis capability 
to model laminated composites for global structural response. 

After each finite element analysis, the CMC module is called to compute constituent 
stresses. These stresses are compared to failure criteria for evaluation of possible fiber, 
matrix, and interface failure in each ply at each node of the ceramic composite structure. 
Once the state of damage at each node is assessed, a damage index is created to record the 
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damage information at each damaged node. The damage index contains the node number, 
the ply number, and the list of damage criteria that have become activated. Fibers, 
matrices, and interfaces are each assigned their own damage index to record the 
corresponding damage information. When a new failure occurs in a fiber, matrix, or 
interface, its damage index is updated. The properties of fibers, matrices and interfaces are 
degraded according to their damage index. For example, if a fiber is damaged due to 
longitudinal tensile failure criteria, its corresponding modulus (£//) is degraded to almost 
zero. 

If there is no damage after a load increment, the structure is considered to be in 
equilibrium and an additional load increment is applied to start the next analysis cycle. At 
each cycle, the process of checking the failure criteria, updating the finite element mesh, 
and updating the material properties is repeated with possible damage growth and 
propagation. Consecutive analysis cycles are continued until structural collapse is 
simulated. 

In order to evaluate the composite fracture toughness, computation of a damage energy 
release rate (DERR) during progressive fracture is needed. DERR is defined as the ratio of 
incremental damage energy to the corresponding incremental damage volume that is 
generated, 

DamageEnergy = ZXM?, /E fj h+EE(°- 5S i , /E «, K /E >, K 
< i i i ‘ i 

Where 

S F j is the fiber strength, j can be 1 1, 22, 33, 12, 23, 13. 
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Efj is the fiber modulus, j can be 1 1, 22, 33, 12, 23, 13. 


Vfj is the volume corresponding to the failed fiber. 

SMj is the matrix strength, j can be 1 1, 22, 33, 12, 23, 13. 

Euj is the matrix modulus, j can be 1 1, 22, 33, 12, 23, 13. 

V Mi is the volume corresponding to the failed matrix. 

Sij is the interface strength, j can be 1 1, 22, 33, 12, 23, 13. 

Sij is the interface modulus, j can be 1 1, 22, 33, 12, 23, 13. 

Vij is the volume corresponding to the failed interface. 

The damage energy release rate is defined as: 

DERR = Damage Energy / Damage Volume 


7.5 Examples and Discussion 

In order to demonstrate the capability to simulate progressive fracture of ceramic matrix 
composites, seven NICALON/1723 CMC specimens tested by Coker and Ashbaugh 
(1992) were simulated. These test specimens were cut from two composite panels, 88C23 
and 89C03, as described by Coker and Ashbaugh (1992). The composite was cross-plied 
with a [(0/90 ) 3 ] s lay-up. The specimen specifications are defined in figure 7.4 and the 
dimensions for the specimen geometry are presented in table 7.1. The fiber volume 
fractions for these specimens depended on the size and location of the sample area, and 
varied from 42 to 48 percent. The average fiber volume fraction of 45% was taken as the 
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fiber volume fraction for simulations in this chapter. The properties for fibers, matrices, 
and interfaces are shown in tables 7.2, 7.3, and 7.4. The thickness of interface is assumed 
to be 2 percent of fiber diameter. 
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B: Thickness 

W: Distance between load line and back face 
H: Half height of specimen 
a: Distance between load line and notch tip 
x: Distance between load line and front face 


Figure 7.4 Definition of Specimen Dimensions 
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Table 7. 1 Specimen Dimensions 


Specimen No. 

B (mm) 

W (mm) 

H (mm) 

X (mm) 

a (mm) 

89C03-2 

3.12 

19.97 

11.96 

5.14 

10.43 

89C03-4 

3.03 

19.96 

11.97 

5.30 

10.39 

89C03-5 

3.08 

20.04 

11.94 

5.09 

9.96 

89C03-7 

3.14 

20.17 

11.95 

4.95 

10.25 

89C03-8 

2.96 

19.81 

12.03 

5.26 

9.91 

89C03-10 

3.10 

40.03 

23.95 

10.04 

19.86 

88C23-7 

2.57 

40.06 

23.92 

10.11 

19.95 


Table 7.2 NICALON fiber properties 


Normal modulus (11) 

2.900E+07 (psi) 

Normal modulus (22) 

2.900E+07 (psi) 

Poisson's ratio (12) 

0.300 

Poisson's ratio (23) 

0.300 

Shear modulus (12) 

1.120E+07 (psi) 

Shear modulus (23) 

1.120E+07 (psi) 

Thermal Expansion Coefficient (11) 

0.210E-05 (in/in/F) 

Thermal Expansion Coefficient (22) 

0.230E-05 (in/in/F) 

Thermal conductivity (11) 

0.208E+01 (BTU/hour/F/in) 

Thermal conductivity (22) 

0.340E+00 (BTU/hour/F/in) 

Heat capacity 

0.246E+00 (BTU/lb) 

Tension strength (11) 

0.350E+06 (psi) 

Compression strength (11) 

0.350E+06 (psi) 

Tension strength (22) 

0.350E+06 (psi) 

Compression strength (22) 

0.350E+06 (psi) 

Torsion strength (12) 

0.240E+06 (psi) 

Torsion strength (23) 

0.240E+06 (psi) 
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Table 7.3 1723 Glass-Matrix properties 


Normal modulus 

0.127E+08 (psi) 

Poisson's ratio 

0.220 

Thermal Expansion Coefficient 

0.240E-05 (in/in/F ) 

Thermal conductivity 

1.167 (BTU/hour/F/in) 

Heat capacity 

0.1563 (BTU/lb) 

Tension strength 

10000.0 (psi) 

Compression strength 

15000.0 (psi) 


Shear strength 8000 (psi) 


Table 7.4 Interface properties 


Normal modulus 

1.590E+06 (psi) 

Poisson's ratio 

0.150 

Thermal Expansion Coefficient 

0.240E-05 (in/in/F) 

Thermal conductivity 

1.20 (BTU/hour/F/in) 

Heat capacity 

0.120 (BTU/lb) 

Tension strength 

0.800E+04 (psi) 

Compression strength 

0. 145E+05 (psi) 

Shear strength 

0.800E+04 (psi) 


A computational model of the specimen contains 439 nodes and 384 elements as shown in 
Figure 7-5. The pin holes were not modeled in the finite element representation of the 
specimen to enable nodal support and loading. The finite element model was configured 
to have a node point at the center of each pin hole. One of the load points (node 371) was 
restrained in all degrees of freedom except for 0 Z . The other load point (node 382) was 
restrained only in D x , D z , 0 X , 9y directions, but allowed motion in Y direction and rotation 
about the Z-axis. A concentrated tensile load was applied at node 382 in the Y direction. 
The load was increased gradually. 

Since all the specimens had similar fracture propagation characteristics, only one 
specimen, 89C03-10, is discussed in detail here. For specimen 89C03-10, the load started 
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from 44N (101b.). When the load was under 262N (59 lb.), there was no damage in the 
specimen. As the load increased to 262N (591b.), the matrix in ply 4 started to fracture at 
the notch tip (node 179) due to transverse tension. With the additional increase of the 
loading, matrix fractures expanded to other plies. When the load increased to 333N 
(751b.), fibers in plies 1,2,4,6,7,9,10 at node 179 were damaged due to longitudinal 
compressive failures. When the load increased to 475N (107 lb.), damage expanded to 
adjacent nodes, 162, 178, 179, 180. When the load increased to 853N (192 lb.), damage 
propagated to nodes 127 - 129, 143 - 147, 160 - 164, 177 - 181. When the load increased 
to 1044N (2351b.), nodes 196 and 197 were completely fractured, and partial-thickness 
damage was propagated to 54 nodes. When the load increased to 147 IN (3311b.), damage 
expanded to 82 nodes, but still only two nodes, 196 and 197, were completely fractured. 
The 331 lb. load was the maximum equilibrium load. After this load, the specimen entered 
an unstable .fracture propagation stage. The load could not be increased above 331 lbs 
without fracturing the specimen. When the load was increased to 1475 N (332 lbs.), the 
specimen broke into two pieces. Figure 7-6 shows the finite element model immediately 
before the structural fracture of this specimen. Figure 7.7 shows damage energy release 
rates with loading. After damage initiation, low points in the DERR indicate reduced 
structural resistance against damage propagation. On the other hand, high points in the 
DERR correspond to structural resistance stages against damage progression. Figure 7.7 
indicates that the specimen demonstrates significant resistance to damage initiation at 
approximately 0.3 kN and also resists damage propagation at a load of 1.3 kN. This 
resistance is confirmed by figure 7.9 where damage progression rate is smaller at this load. 
Figure 7.8 shows the total damage energy release rates. The total DERR is defined as the 
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ratio of the sum of all damage energies to the total damage volume, accumulated from 
damage initiation to the simulated load level. The total DERR presents a smoother overall 
picture of structural fracture resistance characteristics of a specimen. Figure 7.9 shows 
the percent damage volume versus loading. A steep increase in the damage level indicates 
the occurrence of critical damage processes. 

Table 7.5 shows peak loads for both test and simulation. The test results were provided by 
Coker et al (1992). From Table 7.5, it can be seen that simulation results agree with test 
results very well except for specimens 89C03-7 and 89C03-8. From table 7.1, it can be 
seen that the geometry and size of these specimens are close to those of specimens 89C03- 
2, 89C03-4, and 89C03-5. Therefore it seems the peak loads for these specimens should 
also be close as the simulation shows. Test results for specimen 89C03-7 and 89C03-8 
give peak loads that seem too low. This may be explained as the scattering of test data. 
Compared to PMC materials, the properties and peak load of continuous fiber reinforced 
CMC structures are more difficult to predict since CMC materials behave in a brittle 
fashion that increases the uncertainties of the failure load levels. Test results for CMC 
structures usually show more scattering than those of PMC structures. Small variations in 
the fiber and matrix constituent proprieties may cause a significant change in the structural 
fracture load. Therefore, more accurate representations of the composite constituent 
materials are required for CMC micromechanics analysis. 

Figures 7.10 through 7.16 show load-displacement relationships for the seven specimens. 
The displacements shown in these figures are the crack mouth opening displacements 
(CMOD). It should be noted that, simulations are load-controlled, whereas experimental 
results are from displacement-controlled tests. For the load-controlled simulations, the 
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load is increased gradually until the specimen is completely failed. During a load 
controlled simulation there is no unloading process that is taken into account. Therefore, 
after the load reaches its peak value, no additional load can be applied to the specimen and 
the specimen fractures or collapses very quickly. However, during a displacement- 
controlled test, the displacement increases gradually instead of the load. When the 
specimen reaches its maximum strength, the load also reaches its maximum value-or the 
peak load. After the peak load, global structural stiffness degrades significantly. 
Accordingly, the load which can be measured as a reaction under the imposed 
displacement is decreased. The unloading process after the peak load can be plotted in a 
load-displacement relationship for a displacement controlled test. 


Table 7.5 Peak Loads 


Specimen 

Test peak 
loads (KN) 

Simulation peak loads 
(KN) 

Relative Error 
(%) 

89C03-2 

0.799 

0.813 

1.7 

89C03-4 

0.822 

0.800 

2.6 

89C03-5 

0.865 

0.892 

3.1 

89C03-7 

0.631 

0.862 

36.6 

89C03-8 

0.741 

0.866 

16.8 

89C03-10 

1.226 

1.191 

2.8 

88C23-7 

1.457 

1.471 

1.0 
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Total DERR (Mpa) DERR (Mpa) 



Figure 7.7 Damage Energy Release Rate with Loading for Specimen 89C03-10 



Figure 7.8 Total Damage Energy Release Rate with Loading for Specimen 89C03-10 
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Figure 7.9 Damage Propagation with Load for Specimen 89C03-10 



Figure 7.10 Load Displacement Relationship of Specimen 89(303-10 
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Figure 7. 12 Load Displacement Relationship of Specimen 89C03-4 
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Figure 7.13 Load Displacement Relationship of Specimen 89C03-5 



Figure 7. 14 Load Displacement Relationship of Specimen 89(703-7 
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Figure 7.15 Load Displacement Relationship of Specimen 89C03-8 



Figure 7. 16 Load Displacement Relationship of Specimen 89C23-7 
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7.6 Conclusion 


The significant results from this investigation in which computational simulation was used 
to evaluate damage growth and propagation to fracture for SiC/1723 ceramic matrix 
composite specimens are as follows: 

1 . Computational simulation, with the use of established composite mechanics and finite 
element modules, can be used to predict the influence of an existing notch, as well as 
loading, on the progressive fracture of ceramic matrix fiber composite specimens. 

2. Computational simulation adequately tracks the damage growth and subsequent 
propagation to fracture for composite compact tension specimens. 

3. Computational simulation can be used prior to testing to identify locations and modes 
of composite damage that need be monitored by proper instrumentation and inspection of 
the specimen during a laboratory experiment. 

4. Interpretation of experimental data can be facilitated significantly by detailed results 
from a computational simulation. 

5. Computational simulation provides detailed information on damage initiation and 
progression mechanisms, as well as identifying sensitive material parameters affecting 
structural fracture. 

6. The demonstrated procedure is flexible and applicable to all types of constituent 
materials, structural geometry, and loading. Hybrid composite structures composites 
containing homogeneous materials such as metallic layers, as well as binary composites 
can be simulated. 
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7. Fracture toughness parameters such as the structural fracture load are identifiable for 
any specimen or structure by the demonstrated method. 

8. Computational simulation represents a new global approach that may be used for 
progressive damage and fracture assessment in design investigations. 


NASA/CR— 2005-213963 


128 



Chapter 8 

Conclusion 

Predictions of damage growth, and propagation to fracture are needed to evaluate damage 
tolerance, safety, and reliability of composite structures. Damage stability is influenced by 
constituent material properties and global factors, such as structural geometry and 
boundary conditions. The interaction of these factors, further complicated by the 
numerous possibilities of material combinations, composite geometry, fiber orientations, 
and loading conditions, renders the assessment of composite damage progression very 
complex. This complexity makes it very expensive and time consuming to identify and 
isolate all significant parameters from tests. Computerized simulation can reduce 
experimental testing requirements and shorten the design time dramatically. This research 
provided a computational tool for the purpose of progressive fracture and damage 
tolerance assessment for three-dimensionally reinforced fiber composite structures. The 
demonstrated, and verified computer software code was developed by enhancing an 
existing software code, CODSTRAN, with mathematical models and associated software 
codes that provide the capability to predict the effects of mechanical and thermal loads for 
the three dimensional polymer/ceramic matrix composite structures. The tasks carried out 
to accomplish this effort were as follows: 

1. A mathematical model and associated software code were developed to evaluate 
elastic properties and stress limits of three dimensionally reinforced fiber composites. 
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2. A mathematical model and associated software code were developed to 
computationally simulate thermo-mechanical loading of three-dimensional woven 
composite structures. 

3. A mathematical model and associated software code were developed to predict 
damage initiation, propagation, and final failure under mechanical/ thermal loads for 
composite structures. 

4. A mathematical model and associated software code were developed to predict 
material properties (strength, modulus) degradation relative to service conditions. 

5. The developed software simulation results were compared with existing experimental 
test results to validate evaluation of durability and robustness of the computational 
approach for three dimensionally reinforced fiber composites. 

A significant feature of this research is that it takes into account simulation of crack 
propagation in 3D polymer/ceramic matrix materials. This can not be accomplished by the 
use of conventional FEM analysis and laminate theory since classical laminate theory does 
not consider through the thickness aspects of out-of-plane fiber reinforcements in 
composite materials, and FEM analysis does not consider degradation of material 
properties under loading. Further work can be carried out to improve and enhance the 
current research in the following areas: 

1. Progressive damage and fracture of a composite structure subjected to cyclic loading. 

2. Prediction of property degradation for composite structures due tp thermal loading. 

3. Prediction of creep rapture of composite structures. 

4. Probabilistic simulation of fracture to quantify reliability of composite structures. 
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Appendix A: Constituent Material Properties 


The following properties were used in the composite material databank to represent the 
fiber and matrix constituents: 


ASW4 Graphite Fiber Properties: 

Number of fibers per end = 10000 

Fiber diameter = 0.00762 mm (0.300E-3 in) 

Fiber Density = 4.04E-7 Kg/m (0.063 lb/in 3 ) 

Longitudinal normal modulus = 235 GPa (34.08E+6 psi) 

Transverse normal modulus = 17.0 GPa (2.47E+6 psi) 

Poisson's ratio (p. 12 ) = 0.25 
Poisson's ratio (fj.23) =0.27 
Shear modulus (G 12 ) = 55.1 GPa (7.98E+6 psi) 

Shear modulus (G 23 ) = 6.90 GPa (1.00E+6 psi) 

Longitudinal thermal expansion coefficient = -1.0E-6/°C (- 0.55E-6/°F ) 

Transverse thermal expansion coefficient = 1.0E-5/°C (0.56E-5/°F) 

Longitudinal heat conductivity = 0 .302 J-m/hr/m /°C (4.03 BTU -in/hr/in /°F) 

Transverse heat conductivity = 0.0302 J-m/hr/m 2 /°C (0.403 BTU -in/hr/in /°F) 

Heat capacity = 712 J/Kg/°C (0. 1 7 BTU/lb/°F) 

Tensile strength = 3, 723 MPa (540 ksi) 

Compressive strength = 3,351 MPa (486 ksi) 


SGLW Glass Fiber Properties: 

Number of fibers per end = 204 

Fiber diameter = 0.00914 mm (0.360E-3 in) 

Fiber Density = 5.77E-7 Kg/m (0.090 lb/in ) 
Longitudinal normal modulus = 84.8 GPa (12.3E+6 psi) 
Transverse normal modulus = 84.8 GPa (12.3E+6 psi) 
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Poisson’s ratio (p. 12 ) = 0.22 
Poisson's ratio (H 23 ) = 0.22 
Shear modulus (G 12 ) = 34.8 GPa (5.04E+6 psi) 

Shear modulus (G 23 ) = 34.8GPa (5.04E+6 psi) 

Longitudinal thermal expansion coefficient = 0 ,504E-5/°C (0.280E-5/°F) 

Transverse thermal expansion coefficient = 0 ,504E-5/°C (0.280E-5/°F) 
Longitudinal heat conductivity = 3.90E-3 J-m/hr/m 2 /°C (5.208E-2 BTU -in/hr/in /°F) 
Transverse heat conductivity = 3.90E-3 J-m/hr/m 2 /°C (5.208E-2 BTU -in/hr/in 2 /°F) 
Heat capacity = 712 J/Kg/°C (0.17 BTU/lb/°F) 

Tensile strength = 2482 MPa (360 ksi) 

Compressive strength = 2069 MPa (300 ksi) 


Dow Tactix 138 Epoxy resin with H41 hardener Matrix Properties: 

Matrix density = 3.35E-7 Kg/m (0.0450 Ib/in) 

Normal modulus = 2.99 GPa (435 ksi) 

Poisson's ratio = 0.300 

Coefficient of thermal expansion = 0.72E-4/°C (0.40E-4/°F) 

Heat conductivity = 8.681E-3 BTU -in/hr /in /°F 
Heat capacity = 0.25 BTU/lb/°F 

Tensile strength = 85.0 MPa (12.3 ksi) 

Compressive strength = 423 MPa (61.3 ksi) 

Shear strength = 147 MPa (8.17 ksi) 

Allowable tensile strain = 0.02 

Allowable compressive strain = 0.05 

Allowable shear strain = 0 .04 

Allowable torsional strain = 0 .04 

Void conductivity = 16.8 J-m/hr/m 2 /°C (0.225 BTU -in/hr/in /°F) 
Glass transition temperature = 216°C (420°F) 
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Appendix B: Comparison of Composite Properties Computed by the 3-D 
Composite Mechanics (3D-CM) with Cox (1995) and Experimental 
Data 

The following properties are computed by the 3-D Composite Mechanics (3D-CM): 

Exx= Elastic modulus in stuffer direction 

Eyy= Elastic modulus in fill er direction 

Ezz= Elastic modulus in normal direction 

Gyz= Shear modulus in the filler-normal plane 

Gzx= Shear modulus in the stuffer-normal plane 

Gxy= Shear modulus in the stuffer-filler plane 

NUxy= Poisson's ratio in the stuffer-filler plane 

NUyz= Poisson's ratio in the filler-normal plane 

NUxz= Poisson's ratio in the stuffer-normal plane 

SWxxTf= Tensile strength in the stuffer direction based on fiber stress 

SWxxTm= Tensile strength in the stuffer direction based on matrix stress 

SWxxCf= Compressive strength in the stuffer direction based on fiber stress 

SWxxCm= Compressive strength in the stuffer direction based on matrix stress 

SWyyTf= Tensile strength in the filler direction based on fiber stress 

SWyyTm= Tensile strength in the filler direction based on matrix stress 

SWyyCf= Compressive strength in the filler direction based on fiber stress 

SWyyCm= Compressive strength in the filler direction based on matrix stress 

SWzzTf= Tensile strength in the normal direction based on fiber stress 

SWzzTm= Tensile strength in the normal direction based on matrix stress 

SWzzCf= Compressive strength in the normal direction based on fiber stress 

SWzzCm= Compressive strength in the normal direction based on matrix stress 

SWyzS= Shear strength in the filler-normal plane 

SWzxS= Shear strength in the stuffer-normal plane 

SWxyS= Shear strength in the stuffer-filler plane 
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Specimen 1 (LL1) Woven Composite 


Property 

3D-CM 

Cox 

Experiment 

Exx (GPa) 

31.18 

36.8 

30+/-6 

Eyy (GPa) 

32.09 

38.7 


Ezz (GPa) 

8.22 

9.0 

5.7 

Gyz (GPa) 

2.2 

2.1 


Gzx (GPa) 

4.64 

6.0 


Gxy (GPa) 

2.63 

2.3 


Nuxy 

0.0468 

0.023 

0.024 

NUyz 

0.353 

0.216 


NUxz 

0.494 

0.207 

0.22 

SWxxTf(GPa) 

0.631 



SWxxTm(GPa) 

0.0558 



SWxxCf (GPa) 

0.5679 



SWxxCm(GPa) 

0.1954 



SWyyTf (GPa) 

0.5435 



SWyyTm(GPa) 

0.3959 



SWyyCf (GPa) 

0.4892 



SWyyCm(GPa) 

0.1973 



SWzzTf (GPa) 

0.1284 



SWzzTm(GPa) 

0.0821 



SWzzCf(GPa) 

0.1156 



SWzzCih(GPa) 

0.3259 



SWyzS (GPa) 

0.1152 



SWzxS(GPa) 

0.2071 



SWxyS (GPa) 

0.1208 
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Specimen 2 (LL2) Woven Composite 


Property 

3D-CM 

Cox 

Experiment 

Exx (GPa) 

32.79 

34.9 

28.5 

Eyy (GPa) 

42.05 

47.6 


Ezz (GPa) 

8.06 

7.0 

5.9 

Gyz (GPa) 

2.44 

2.2 


Gzx (GPa) 

3.10 

3.2 


Gxy (GPa) 

2.89 

2.4 


NUxy 

0.040 

0.027 

0.11 

NUyz 

0.376 

0.310 


NUxz 

0.460 

0.457 

0.50 

SWxxT(GPa) 

0.5833 



SWxxT(GPa) 

0.0513 



SWxxCf (GPa) 

0.5208 



SWxxCm (GPa) 

0.2078 



SWyyTf (GPa) 

0.7446 



SWyyTm (GPa) 

0.3297 



SWyyCf (GPa) 

0.6701 



SWyyCm (GPa) 

0.1643 



SWzzTf (GPa) 

0.0620 



SWzzTm (GPa) 

0.0759 



SWzzCf (GPa) 

0.0517 



SWzzCm (GPa) 

0.3302 



SWyzS (GPa) 

0.1102 



SWzxS (GPa) 

0.1473 



SWxyS (GPa) 

0.1179 
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Specimen 3(LT1) Woven Composite 


Property 

3D-CM 

Cox 

Experiment 

Exx (GPa) 

42.21 

47.3 

27 

Eyy (GPa) 

52.98 

59.5 


Ezz (GPa) 

9.27 

9.4 

8.0 

Gyz (GPa) 

2.78 

2.7 


Gzx (GPa) 

4.53 

5.6 


Gxy (GPa) 

3.38 

3.0 


NUxy 

0.0325 

0.02 

0.048 

NUyz 

0.342 

0.243 


NUxz 

0.470 

0.541 

0.375 

SWxxTf (GPa) 

0.7608 



SWxxTm (GPa) 

0.0511 



SWxxCf (GPa) 

0.6027 



SWxxCm (GPa) 

0.2069 



SWyyTf (GPa) 

0.8744 



SWyyTm (GPa) 

0.0345 



SWyyCf (GPa) 

0.6927 



SWyyCm (GPa) 

0.1722 



SWzzTf (GPa) 

0.0998 



SWzzTm (GPa) 

0.0777 



SWzzCf (GPa) 

0.0790 



SWzzCm (GPa) 

0.3392 



SWyzS (GPa) 

0.1100 



SWzxS (GPa) 

0.1751 



SWxyS (GPa) 

0.1208 
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Specimen 4(LT2) Woven Composite 


Property 

3D-CM 

Cox 

Experiment 

Exx (GPa) 

42.68 

43.5 

39 

Eyy (GPa) 

50.56 

51.6 


Ezz (GPa) 

8.54 

7.0 

7.9 

Gyz (GPa) 

2.69 

2.4 


Gzx (GPa) 

3.15 

3.1 


Gxy (GPa) 

3.31 

2.6 


NUxy 

0.036 

0.027 

0.021 

NUyz 

0.379 

0.325 


NUxz 

0.436 

0.37 

0.37 

SWxxTf (GPa) 

0.7377 



SWxxTm (GPa) 

0.0466 



SWxxCf (GPa) 

0.6012 



SWxxCm (GPa) 

0.1969 



SWyyTf (GPa) 

0.8363 



SWyyTm (GPa) 

0.0346 



SWyyCf (GPa) 

0.6805 



SWyyCm (GPa) 

0.1727 



SWzzTf (GPa) 

0.0544 



SWzzTm (GPa) 

0.0748 



SWzzCf (GPa) 

0.0453 



SWzzCm (GPa) 

0.3374 



SWyzS (GPa) 

0.1083 



SWzxS (GPa) 

0.1426 



SWxyS (GPa) 

0.1187 
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Specimen 5(L-0) Woven Composite 


Property 

3D-CM 

Cox 

Experiment 

Exx (GPa) 

47.34 

51.9 

30+/-2 

Eyy (GPa) 

57.50 

63.9 

45.5+/- 1.5 

Ezz (GPa) 

13.97 

13.7 

7.0+/- 1 

Gyz (GPa) 

2.91 

2.7 


Gzx (GPa) 

2.73 

2.8 


Gxy (GPa) 

3.51 

3.1 


NUxy 

0.0446 

0.034 

0.053 

NUyz 

0.236 

0.183 


NUxz 

0.242 

0.184 

0.49 

SWxxTf (GPa) 

0.7499 



SWxxTm (GPa) 

0.0409 



SWxxCf (GPa) 

0.5759 



SWxxCm (GPa) 

0.2037 



SWyyTf (GPa) 

0.9423 



SWyyTm (GPa) 

0.0334 



SWyyCf (GPa) 

0.7238 



SWyyCm (GPa) 

0.1664 



SWzzTf (GPa) 

0.1079 



SWzzTm (GPa) 

0.0658 



SWzzCf (GPa) 

0.0828 



SWzzCm (GPa) 

0.3282 



SWyzS (GPa) 

0.1048 



SWzxS (GPa) 

0.1040 



SWxyS (GPa) 

0.1162 
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Specimen 6(HL1) Woven Composite 


Property 

3D-CM 

Cox 

Experiment 

Exx (GPa) 

87.31 

91.5 

85+/-8 

Eyy (GPa) 

54.34 

56.2 

43.8 

Ezz (GPa) 

12.58 

12.1 

16+/-2 

Gyz (GPa) 

4.47 

4.1 


Gzx (GPa) 

6.75 

7.1 


Gxy (GPa) 

5.99 

5.4 

6.2 

NUxy 

0.0437 

0.034 

0.061 

NUyz 

0.324 

0.286 


NUxz 

0.422 

0.456 


SWxxTf (GPa) 

1.1430 



SWxxTm (GPa) 

0.0357 



SWxxCf (GPa) 

0.9021 



SWxxCm (GPa) 

0.1462 



SWyyTf (GPa) 

0.7848 



SWyyTm (GPa) 

0.0486 



SWyyCf(GPa) 

0.4919 



SWyyCm (GPa) 

0.2424 



SWzzTf (GPa) 

0.0842 



SWzzTm (GPa) 

0.0789 



SWzzCf (GPa) 

0.0528 



SWzzCm (GPa) 

0.3618 



SWyzS (GPa) 

0.1130 



SWzxS (GPa) 

0.1565 



SWxyS (GPa) 

0.1264 
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Specimen 7(HL2) Woven Composite 


Property 

3D-CM 

Cox 

Experiment 

Exx (GPa) 

75.32 

81.2 

80 

Eyy (GPa) 

51.69 

55 

42.3 

Ezz (GPa) 

11.21 

10.2 

14.0 

Gyz (GPa) 

3.61 

3.6 


Gzx (GPa) 

5.26 

5.3 


Gxy (GPa) 

5.16 

4.6 

5.8 

NUxy 

0.0439 

0.035 

0.13 

NUyz 

0.354 

0.298 


NUxz 

0.420 

0.425 

0.45+/-0.05 

SWxxTf (GPa) 

1.2495 



SWxxTm (GPa) 

0.0335 



SWxxCf (GPa) 

0.8527 



SWxxCm (GPa) 

0.1479 



SWyyTf (GPa) 

0.7776 



SWyyTm (GPa) 

0.0451 



SWyyCf (GPa) 

0.5307 



SWyyCrn (GPa) 

0.2250 



SWzzTf (GPa) 

0.0466 



SWzzTm (GPa) 

0.0754 



SWzzCf (GPa) 

0.0318 



SWzzCm (GPa) 

0.3567 



SWyzS (GPa) 

0.1102 



SWzxS (GPa) 

0.1364 



SWxyS (GPa) 

0.1221 
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Specimen 8(HT1) Woven Composite 


Property 

3D-CM 

Cox 

Experiment 

Exx (GPa) 

84.24 

88.6 

79 

Eyy (GPa) 

52.48 

54.4 

42.5 

Ezz (GPa) 

12.80 

12.8 

13.8 

Gyz (GPa) 

3.99 

4.0 


Gzx (GPa) 

7.12 

7.8 


Gxy (GPa) 

5.84 

5.3 

5.6 

NUxy 

0.0444 

0.033 

0.054 

NUyz 

0.316 

0.248 


NUxz 

0.430 

0.486 


SWxxTf (GPa) 

1.415 



SWxxTm (GPa) 

0.0386 



SWxxCf (GPa) 

0.8948 



SWxxCm (GPa) 

0.1499- 



SWyyTf (GPa) 

0.7554 



SWyyTm (GPa) 

0.0492 



SWyyCf (GPa) 

0.4777 



SWyyCm (GPa) 

0.2451 



SWzzTf (GPa) 

0.1118 



SWzzTm (GPa) 

0.0806 



SWzzCf (GPa) 

0.0707 



SWzzCm (GPa) 

0.3591 



SWyzS (GPa) 

0.1142 



SWzxS (GPa) 

0.1724 



SWxyS (GPa) 

0.1271 
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Specimen 9(HT2) Woven Composite 


Property 

3D-CM 

Cox 

Experiment 

Exx (GPa) 

80.26 

85.1 

72 

Eyy (GPa) 

54.91 

57.6 

45.8 

Ezz (GPa) 

11.92 

11.2 

13.9 

Gyz (GPa) 

3.88 

3.9 


Gzx (GPa) 

5.98 

6.2 


Gxy (GPa) 

5.59 

5.0 

5.7 

NUxy 

0.0423 

0.033 

0.097 

NUyz 

0.337 

0.280 


NUxz 

0.422 

0.443 


SWxxTf (GPa) 

1.323 



SWxxTm (GPa) 

0.0355 



SWxxCf (GPa) 

0.8602 



SWxxCm (GPa) 

0,1514 



SWyyTf (GPa) 

0.8133 



SWyyTm (GPa) 

0.0460 



SWyyCf (GPa) 

0.5285 



SWyyCm (GPa) 

0.2293 



SWzzTf (GPa) 

0.0650 



SWzzTm (GPa) 

0.0772 



SWzzCf (GPa) 

0.0422 



SWzzCm (GPa) 

0.3592 



SWyzS (GPa) 

0.1114 



SWzxS (GPa) 

0.1462 



SWxyS (GPa) 

0.1243 
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Specimen lO(HOl) Woven Composite 


Property 

3D-CM 

Cox 

Experiment 

Exx (GPa) 

89.08 

93.1 

88 

Eyy (GPa) 

53.71 

56.4 

39.9 

Ezz (GPa) 

16.14 

17.3 

15.4 

Gyz (GPa) 

3.41 

4.1 


Gzx (GPa) 

4.12 

4.7 


Gxy (GPa) 

4.89 

5.4 

5.0 

NUxy 

0.0539 

0.051 

0.055 

NUyz 

0.230 

0.192 


NUxz 

0.221 

0.190 


SWxxTf (GPa) 

1.405 



SWxxTm (GPa) 

0.0284 



SWxxCf (GPa) 

0.8815 



SWxxCm (GPa) 

0.1416 



SWyyTf (GPa) 

0.7835 



SWyyTm (GPa) 

0.0481 



SWyyCf (GPa) 

0.4913 



SWyyCm (GPa) 

0.2400 



SWzzTf (GPa) 

0.1129 



SWzzTm (GPa) 

0.0694 



SWzzCf (GPa) 

0.0708 



SWzzCm (GPa) 

0.3460 



SWyzS (GPa) 

0.1058 



SWzxS (GPa) 

0.1049 



SWxyS (GPa) 

0.1222 
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Specimen 1 1(H02) Woven Composite 


Property 

3D-CM 

Cox 

Experiment 

Exx (GPa) 

80.54 

83.8 

69+/-5 

Eyy (GPa) 

52.68 

55.9 

41.6 

Ezz (GPa) 

17.97 

20.4 

22.3 

Gyz (GPa) 

3.30 

4.0 


Gzx (GPa) 

3.84 

4.4 


Gxy (GPa) 

4.52 

4.9 


NUxy 

0.0559 

0.052 

0.07 

NUyz 

0.203 

0.158 


NUxz 

0.196 

0.157 


SWxxTf (GPa) 

1.2761 



SWxxTm (GPa) 

0.0305 



SWxxCf (GPa) 

0.8275 



SWxxCm (GPa) 

0.1523 



SWyyTf (GPa) 

0.7793 



SWyyTm (GPa) 

0.0468 



SWyyCf (GPa) 

0.5053 



SWyyCm (GPa) 

0.2335 



SWzzTf (GPa) 

0.1523 



SWzzTm (GPa) 

0.0674 



SWzzCf (GPa) 

0.0988 



SWzzCm (GPa) 

0.3360 



SWyzS (GPa) 

0.1057 



SWzxS (GPa) 

0.1045 



SWxyS (GPa) 

0 . 1207 
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degradation, from damage initiation to collapse of structures, in the 3D woven structures are simulated for the first time. 
Three dimensional woven composite specimens with various woven patterns under different loading conditions, such as 
tension, compression, bending, and shear are simulated in the validation process of this research. Damage initiation, 
growth, accumulation, and propagation to fracture are included in these simulations. 
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